list.of.packages <- c("stargazer","gdata","spdep","maptools", "GISTools","geosphere","rgdal","rgeos","AER","spdep","raster","xtable","mvtnorm","stringr","cem","MatchIt","RCurl","mgcv")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]; if(length(new.packages)){install.packages(new.packages)}
lapply(list.of.packages, require, character.only = TRUE)

rm(list=ls())
Sys.setlocale('LC_ALL', 'russian')
setwd("~/MyDirectory/ReplicationFiles")
setwd("~/Dropbox/BELARUS/WP/ReplicationFiles/")
source("functions.R")


######################################
# FIGURE 1
######################################

# Load & crop data
rail <- readShapeLines("Shapes/BELARUS_RAIL")
oblasts <- readShapePoly("Shapes/BELARUS_OBLASTS")
load("Shapes/BELARUS_RAYONS_v8.RData")
load("Shapes/BELARUS_VILLAGES_v2.RData")
rail <- crop(as(rail,"SpatialLines"),oblasts)
borders <- readShapePoly("Shapes/April_30_1941")
proj4string(borders) <- CRS("+proj=aea +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs")
proj4string(map) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
borders <- spTransform(borders,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
borders <- gBuffer(borders, byid=TRUE, width=0)
bb <- bbox(map)
bb[1,1] <- bb[1,1]-(bb[1,2]-bb[1,1])/8
bb[1,2] <- bb[1,2]+(bb[1,2]-bb[1,1])/8
bb[2,1] <- bb[2,1]-(bb[2,2]-bb[2,1])/8
bb[2,2] <- bb[2,2]+(bb[2,2]-bb[2,1])/8
WKTs <- paste("POLYGON((",bb[1,1],bb[2,1],",",bb[1,2],bb[2,1],",",bb[1,2],bb[2,2],",",bb[1,1],bb[2,2],",",bb[1,1],bb[2,1],"))")
test <- readWKT(WKTs)
r <- borders
proj4string(test) <- CRS(proj4string(map))
proj4string(r) <- CRS(proj4string(map))
borders2 <- crop(r, test)


# scale points
pnts.size <- function(x,scale=.2){
  ifelse(x<=1,scale*1,
         ifelse(x>1&x<=10,scale*2,
                ifelse(x>10&x<=50,scale*3,
                       ifelse(x>10&x<=100,scale*4,
                              ifelse(x>100&x<=500,scale*5,scale*6)))))}
cexs <- pnts.size(sp.villages$HOUSES_DESTROYED,scale=.15)

# plot: derailment
png("Output/FIG1a.png",width=4,height=3,res=300,units="in")
par(mar=c(0,0,0,0))
plot(test,col="lightblue",border=NA)
plot(borders2,add=T,col="beige",border="grey")
plot(map,add=T,col="beige")
scl <- map$CONVOYS_PROP
#plot(map,col=rgb(1-scl,1-scl,1-scl,.5),add=T,border="lightgray",lwd=.1)
plot(map,col=rgb(scl/2,scl/3,scl,.75),add=T,border="lightgray",lwd=.1)
lines(rail,lty=1,col="black",lwd=3)
lines(rail,lty=1,col="yellow",lwd=2)
lls <- c(bbox(map)[1,1]-2.5,mean(c(bbox(map)[2,1]+1.2,1.1*(bbox(map)[2,2]-.7))))
text(x=lls[1]+2.2,y=lls[2]-.1,labels=c("Partisan violence:"),cex=.7,pos=4)
legend(x=lls[1]+2.3,y=lls[2]-.1,fill=c(rgb(.5,.33,1,.75),rgb(0,0,0,.75),rgb(.25,.167,.5,.75)),legend=c("100% interdiction","100% coercion","50%-50%"),bty="n",cex=.65,ncol=3)
text(bb[1,1]-.1,bb[2,1]+.5,labels="General Government\n(Poland)",pos=4,cex=.5,col="grey")
text(bb[1,1]+4,bb[2,1]+.8,labels="Soviet Union",pos=4,cex=.5,col="grey")
text(bb[1,1]-.1,bb[2,2]-3,labels="Germany",pos=4,cex=.5,col="grey")
legend(x=bb[1,2]-3.5,y=bb[2,2]+.1,lwd=3,col="black",legend="Railroad",bty="n",cex=.7)
legend(x=bb[1,2]-3.5,y=bb[2,2]+.1,lwd=1,col="yellow",legend="Railroad",bty="n",cex=.7)
dev.off()


# plot: reprisals
png("Output/FIG1b.png",width=4,height=3,res=300,units="in")
par(mar=c(0,0,0,0))
plot(test,col="lightblue",border=NA)
plot(borders2,add=T,col="beige",border="grey")
plot(map,add=T,col="beige")
points(sp.villages,pch=16,col=rgb(1,0,0,alpha=.2),cex=cexs)
plot(map,col=NA,border="black",add=T)
text(bb[1,1]-.1,bb[2,1]+.5,labels="General Government\n(Poland)",pos=4,cex=.5,col="grey")
text(bb[1,1]+4,bb[2,1]+.8,labels="Soviet Union",pos=4,cex=.5,col="grey")
text(bb[1,1]-.1,bb[2,2]-3,labels="Germany",pos=4,cex=.5,col="grey")
lls <- c(bbox(map)[1,1]-2.5,mean(c(bbox(map)[2,1]+1.2,1.1*(bbox(map)[2,2]-.7))))
text(x=lls[1]+2.2,y=lls[2]-.1,labels=c("Houses destroyed:"),cex=.7,pos=4)
for(i in 1:6){points(x=lls[1]+6+.85*i,y=lls[2]-.25,cex=pnts.size(c(1,4,19,101,501)[i],scale=.15),pch=16,col="red"); text(x=lls[1]+6+.85*i,y=lls[2]-.05,label=c("1","<10","<100","<500",">500")[i],cex=.5)	}
legend(x=bb[1,2]-3,y=bb[2,2]+.1,fill="beige",border="black",legend="District",bty="n",cex=.7)
dev.off()




######################################
# FIGURE 2
######################################

load("Data/BELARUS_RAYONS_MONTHS_Full.RData")
mat <- rayons.panel
ts.vec <- aggregate(mat[,c("CONVOYS","COERCION","HOUSES_DESTROYED","PEOPLE_KILLED","REPRISAL")],by=list(YRMO=mat$YRMO),FUN=function(x){sum(x,na.rm=T)})
ts.vec$TID <- 1:nrow(ts.vec)

# derail
png("Output/FIG2a.png",width=8,height=2,res=300,units="in")
par(mar=c(2,2,0.5,2))
plot(0,0,col=NA,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",xlim=range(ts.vec$TID),ylim=c(0,500))
axis(1,at=seq(2,max(ts.vec$TID),by=6),labels=paste0(substr(ts.vec$YRMO,5,6),"/",substr(ts.vec$YRMO,1,4))[seq(2,max(ts.vec$TID),by=6)],cex.axis=.75,lwd=2)
axis(1,at=1:max(ts.vec$TID),labels=FALSE,cex.axis=.5)
axis(2,at=seq(0,500,by=100),labels=seq(0,500,by=100),cex.axis=.75,lwd=2,col="purple",col.axis="purple")
axis(2,at=seq(0,500,by=50),labels=FALSE,col="purple")
segments(x0=ts.vec$TID-.2,x1=ts.vec$TID-.2,y0=0,y1=ts.vec$CONVOYS,col="purple",lwd=2)
par(new=TRUE)
plot(0,0,col=NA,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",xlim=range(ts.vec$TID),ylim=c(0,40))
segments(x0=ts.vec$TID+.2,x1=ts.vec$TID+.2,y0=0,y1=ts.vec$COERCION,col="black",lwd=2)
axis(4,at=seq(0,40,by=10),labels=seq(0,40,by=10),cex.axis=.75,lwd=2,col="black")
axis(4,at=seq(0,40,by=5),labels=FALSE)
legend("topleft",fill=c("purple"),legend=c("Interdiction (convoys derailed)"),bty="n")
legend("topright",fill=c("black"),legend=c("Coercion (garrison attacks)"),bty="n")
dev.off()

# reprisals
png("Output/FIG2b.png",width=8,height=2,res=300,units="in")
par(mar=c(2,2,0.5,2))
plot(0,0,col=NA,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",xlim=range(ts.vec$TID),ylim=c(0,50000))
axis(1,at=seq(2,max(ts.vec$TID),by=6),labels=paste0(substr(ts.vec$YRMO,5,6),"/",substr(ts.vec$YRMO,1,4))[seq(2,max(ts.vec$TID),by=6)],cex.axis=.75,lwd=2)
axis(1,at=1:max(ts.vec$TID),labels=FALSE,cex.axis=.5)
axis(2,at=seq(0,50000,by=10000),labels=seq(0,50000,by=10000),cex.axis=.75,lwd=2)
axis(2,at=seq(0,50000,by=2000),labels=FALSE)
segments(x0=ts.vec$TID-.2,x1=ts.vec$TID-.2,y0=0,y1=ts.vec$HOUSES_DESTROYED,col="red",lwd=2)
segments(x0=ts.vec$TID+.2,x1=ts.vec$TID+.2,y0=0,y1=ts.vec$PEOPLE_KILLED,col="blue",lwd=2)
legend("topleft",fill=c("red","blue"),legend=c("Houses destroyed","People killed"),bty="n")
dev.off()



######################################
# TABLE 1
######################################


load("Data/BELARUS_RAYONS_MONTHS_Full.RData")
mat <- rayons.panel
varz <- c("PEOPLE_KILLED","HOUSES_DESTROYED","REPRISAL")

crosstab2 <- lapply(1:3,function(v){
  vard <- paste0(varz[v],"_p1")
  vard_t <- paste0(varz[v],"_t1")
  
  # Post
  y00 <- sum(mat[,vard][mat$CONVOYS==0&mat$COERCION==0],na.rm=T)/sum(mat[,vard],na.rm=T)
  y00p <- mean(mat[,vard][mat$CONVOYS==0&mat$COERCION==0],na.rm=T)
  y10 <- sum(mat[,vard][mat$CONVOYS>0&mat$COERCION==0],na.rm=T)/sum(mat[,vard],na.rm=T)
  y10p <- mean(mat[,vard][mat$CONVOYS>0&mat$COERCION==0],na.rm=T)
  y01 <- sum(mat[,vard][mat$CONVOYS==0&mat$COERCION>0],na.rm=T)/sum(mat[,vard],na.rm=T)
  y01p <- mean(mat[,vard][mat$CONVOYS==0&mat$COERCION>0],na.rm=T)
  y11 <- sum(mat[,vard][mat$CONVOYS>0&mat$COERCION>0],na.rm=T)/sum(mat[,vard],na.rm=T)
  y11p <- mean(mat[,vard][mat$CONVOYS>0&mat$COERCION>0],na.rm=T)
  # Pre
  y00_t <- sum(mat[,vard_t][mat$CONVOYS==0&mat$COERCION==0],na.rm=T)/sum(mat[,vard_t],na.rm=T)
  y00p_t <- mean(mat[,vard_t][mat$CONVOYS==0&mat$COERCION==0],na.rm=T)
  y10_t <- sum(mat[,vard_t][mat$CONVOYS>0&mat$COERCION==0],na.rm=T)/sum(mat[,vard_t],na.rm=T)
  y10p_t <- mean(mat[,vard_t][mat$CONVOYS>0&mat$COERCION==0],na.rm=T)
  y01_t <- sum(mat[,vard_t][mat$CONVOYS==0&mat$COERCION>0],na.rm=T)/sum(mat[,vard_t],na.rm=T)
  y01p_t <- mean(mat[,vard_t][mat$CONVOYS==0&mat$COERCION>0],na.rm=T)
  y11_t <- sum(mat[,vard_t][mat$CONVOYS>0&mat$COERCION>0],na.rm=T)/sum(mat[,vard_t],na.rm=T)
  y11p_t <- mean(mat[,vard_t][mat$CONVOYS>0&mat$COERCION>0],na.rm=T)
  
  # Diff in diff
  mat$DIFF_TEMP <- mat[,vard]-mat[,vard_t]
  vardd <- "DIFF_TEMP"
  y00d <- mean(mat[,vardd][mat$CONVOYS==0&mat$COERCION==0],na.rm=T)
  y10d <- mean(mat[,vardd][mat$CONVOYS>0&mat$COERCION==0],na.rm=T)
  y01d <- mean(mat[,vardd][mat$CONVOYS==0&mat$COERCION>0],na.rm=T)
  y11d <- mean(mat[,vardd][mat$CONVOYS>0&mat$COERCION>0],na.rm=T)
  
  # cross-tab
  crosstab <- data.frame(
    Interdiction0=c(
      paste0(round(y00p,2)),
      paste0(round(y00p_t,2)),
      paste0("(",ifelse(y00d>=0,"+",""),round(y00d,2),")"),
      paste0(round(y01p,2)),
      paste0(round(y01p_t,2)),
      paste0("(",ifelse(y01d>=0,"+",""),round(y01d,2),")")
    ),
    Interdiction1=c(
      paste0(round(y10p,2)),
      paste0(round(y10p_t,2)),
      paste0("(",ifelse(y10d>=0,"+",""),round(y10d,2),")"),
      paste0(round(y11p,2)),
      paste0(round(y11p_t,2)),
      paste0("(",ifelse(y11d>=0,"+",""),round(y11d,2),")")
    ))
  names(crosstab) <- paste0(varz[v],"_",0:1)
  crosstab
  
})
crosstab2 <- do.call(cbind,crosstab2)
xtable(crosstab2)
print.xtable(xtable(crosstab2),type="latex",file=paste0("Output/TABLE1.tex"),include.rownames =TRUE)


######################################
# TABLE 2
######################################


# Village-level
load("Data/BELARUS_VILLAGES_MONTHS_Full.RData")
head(villages.panel)
temp <- data.frame(VID=villages.panel$VID,LONG=villages.panel$LONG,LAT=villages.panel$LAT,DIST2BERLIN=NA)
temp <- temp[!duplicated(temp$VID),]
berlin <- c(13.383333,52.516667)
coords <- temp[,c("LONG","LAT")]
for(i in 1:nrow(coords)){temp[i,"DIST2BERLIN"] <- dist(rbind(berlin,coords[i,]))}
temp$LONG <- temp$LAT <- NULL
X <- villages.panel
X <- merge(X,temp,by="VID",all.x=T,all.y=F)
deg2km <- function(x){x*111.1949}
X$DIST2BERLIN <- deg2km(X$DIST2BERLIN)
varz <- c("PEOPLE_KILLED","HOUSES_DESTROYED","REPRISAL","CONVOYS","COERCION","CONVOYS_PROP","CENT_BTWN_R","DIST_RAIL","URBAN100","OPEN_TERRAIN","ETH_JUD","CONTROL_AREA","DISTLIB_MIN","PEOPLE_PRE","HOUSES_PRE","DIST2BERLIN")
labz <- c("People killed","Houses destroyed","Number of reprisals","Partisan interdiction","Partisan coercion","Partisan interdiction (proportion)","Rail network centrality","Distance to rail station","Urbanization","Open terrain","Percent Jewish","Partisan control","Distance from frontline","People (pre-war)","Houses (pre-war)","Distance to Berlin")
print.xtable(xtable(sumstat(data=X,varz=varz,labs=labz)),file="Output/TABLE2_vil.tex")


# Rayon-level
load("Data/BELARUS_RAYONS_MONTHS_Full.RData")
load("Data/BELARUS_RAYONS_Full_sp.RData")
map_crd <- data.frame(RAYON_ID=map$RAYON_ID,LONG=coordinates(map)[,1],LAT=coordinates(map)[,2])
rayons.panel <- merge(rayons.panel,map_crd,by="RAYON_ID")
temp <- data.frame(RAYON_ID=map$RAYON_ID,DIST2BERLIN=NA)
berlin <- c(13.383333,52.516667)
coords <- coordinates(map)
for(i in 1:nrow(coords)){temp[i,"DIST2BERLIN"] <- dist(rbind(berlin,coords[i,]))}
deg2km <- function(x){x*111.1949}
X <- rayons.panel
X <- merge(X,temp,by="RAYON_ID",all.x=T,all.y=F)
X$DIST2BERLIN <- deg2km(X$DIST2BERLIN)
varz <- c("PEOPLE_KILLED","HOUSES_DESTROYED","REPRISAL","CONVOYS","COERCION","CONVOYS_PROP","CENT_BTWN_R","DIST_RAIL","URBAN100","OPEN_TERRAIN","ETH_JUD","CONTROL_AREA","DISTLIB_MIN","PEOPLE_PRE","HOUSES_PRE","DIST2BERLIN")
labz <- c("People killed","Houses destroyed","Number of reprisals","Partisan interdiction","Partisan coercion","Partisan interdiction (proportion)","Rail network centrality","Distance to rail station","Urbanization","Open terrain","Percent Jewish","Partisan control","Distance from frontline","People (pre-war)","Houses (pre-war)","Distance to Berlin")
print.xtable(xtable(sumstat(data=X,varz=varz,labs=labz)),file="Output/TABLE2.tex")


######################################
# FIGURE 3a
######################################

# load data
load("Data/BELARUS_RAYONS_MONTHS_Full.RData")
# no. of rayons
length(unique(rayons.panel$RAYON_ID))
# no. of time periods
length(unique(rayons.panel$YRMO))
# date range
range(rayons.panel$YRMO)
# N
dim(rayons.panel)

# spatial data
load("Data/BELARUS_RAYONS_Full_sp.RData")
map_crd <- data.frame(RAYON_ID=map$RAYON_ID,LONG=coordinates(map)[,1],LAT=coordinates(map)[,2])
rayons.panel <- merge(rayons.panel,map_crd,by="RAYON_ID")

mod <- gam(PEOPLE_KILLED_p1~as.factor(YRMO)+CONVOYS + COERCION+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN +PEOPLE_PRE+PEOPLE_KILLED+s(LONG,LAT),data=rayons.panel,family="quasipoisson");summary(mod);#qaic(mod)
mod1 <- mod

mod <- gam(PEOPLE_KILLED_p1~as.factor(YRMO)+CONVOYS_PROP +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+PEOPLE_PRE+PEOPLE_KILLED+s(LONG,LAT),data=rayons.panel,family="quasipoisson");summary(mod);#qaic(mod)
mod2 <- mod

mod <- gam(HOUSES_DESTROYED_p1~as.factor(YRMO)+CONVOYS + COERCION +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+HOUSES_PRE+HOUSES_DESTROYED+s(LONG,LAT),data=rayons.panel,family="quasipoisson");summary(mod);#qaic(mod)
mod3 <- mod

mod <- gam(HOUSES_DESTROYED_p1~as.factor(YRMO)+CONVOYS_PROP+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+HOUSES_PRE +HOUSES_DESTROYED+s(LONG,LAT),data=rayons.panel,family="quasipoisson");summary(mod);#qaic(mod)
mod4 <- mod

mod <- gam(REPRISAL_p1~as.factor(YRMO)+CONVOYS +COERCION+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+ N_VILLAGES +REPRISAL+s(LONG,LAT),data=rayons.panel,family="quasipoisson");summary(mod);#qaic(mod)
mod5 <- mod

mod <- gam(REPRISAL_p1~as.factor(YRMO)+CONVOYS_PROP +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+N_VILLAGES+REPRISAL+s(LONG,LAT),data=rayons.panel,family="quasipoisson");summary(mod);#qaic(mod)
mod6 <- mod



# CF plots
png("Output/FIG3.png",height=3,width=12,units="in",res=300)
par(mfrow=c(1,3),mar=c(4,4,.5,.5))
mod <- mod3
contplot_GAM(mod,var="CONVOYS",xlab="Partisan interdiction (convoys derailed)",ylab="Houses destroyed by Germans",cax=1.4,fax=list(YRMO=194310, RAYON_ID=118),col.ci=rgb(.5,.33,1,.75))
contplot_GAM(mod,var="COERCION",xlab="Partisan coercion (garrisons attacked)",ylab="Houses destroyed by Germans",cax=1.4,fax=list(YRMO=194310, RAYON_ID=118))
mod <- mod4
contplot_GAM(mod,var="CONVOYS_PROP",xlab="Partisan interdiction (proportion)",ylab="Houses destroyed by Germans",cax=1.4,fax=list(YRMO=194310, RAYON_ID=118),col.ci=rgb(.5,.33,1,.75))
dev.off()


## CFs
# People killed
preds <- predmanGAM(mod=mod3,var="CONVOYS",vals=c(0,10),quadratic=F,fax=list(YRMO=194310))
paste0(round(preds$y[2]-preds$y[1],2)," (",round(preds$lo[2]-preds$lo[1],2),",",round(preds$up[2]-preds$up[1],2),")")
paste0(round(preds$rr,2)," (",round(preds$rr.lo,2),",",round(preds$rr.up,2),")")

preds <- predmanGAM(mod=mod3,var="COERCION",vals=c(0,1),quadratic=F,fax=list(YRMO=194310))
paste0(round(preds$y[2]-preds$y[1],2)," (",round(preds$lo[2]-preds$lo[1],2),",",round(preds$up[2]-preds$up[1],2),")")
paste0(round(preds$rr,2)," (",round(preds$rr.lo,2),",",round(preds$rr.up,2),")")

preds <- predmanGAM(mod=mod6,var="CONVOYS_PROP",vals=c(0,1),quadratic=F,fax=list(YRMO=194310))
paste0(round(preds$y[2]-preds$y[1],2)," (",round(preds$lo[2]-preds$lo[1],2),",",round(preds$up[2]-preds$up[1],2),")")
paste0(round(preds$rr,2)," (",round(preds$rr.lo,2),",",round(preds$rr.up,2),")")

preds <- predmanGAM(mod=mod1,var="CONTROL_AREA",vals=c(0,.5,1),quadratic=F,fax=list(YRMO=194310))
paste0(round(preds$y[1],2)," (",round(preds$lo[1],2),",",round(preds$up[1],2),")")
paste0(round(preds$y[3],2)," (",round(preds$lo[3],2),",",round(preds$up[3],2),")")
paste0(round(preds$y[2],2)," (",round(preds$lo[2],2),",",round(preds$up[2],2),")")
paste0(round(preds$y[2]-preds$y[1],2)," (",round(preds$lo[2]-preds$lo[1],2),",",round(preds$up[2]-preds$up[1],2),")")
paste0(round(preds$rr,2)," (",round(preds$rr.lo,2),",",round(preds$rr.up,2),")")

preds <- predmanGAM(mod=mod1,var="ETH_JUD",vals=c(1,30),quadratic=F,fax=list(YRMO=194310))
paste0(round(preds$y[2]-preds$y[1],2)," (",round(preds$lo[2]-preds$lo[1],2),",",round(preds$up[2]-preds$up[1],2),")")
paste0(round(preds$rr,2)," (",round(preds$rr.lo,2),",",round(preds$rr.up,2),")")

summary(rayons.panel$DISTLIB_MIN)
preds <- predmanGAM(mod=mod1,var="DISTLIB_MIN",vals=c(77,578),quadratic=F,fax=list(YRMO=194310))
paste0(round(preds$y[2]-preds$y[1],2)," (",round(preds$lo[2]-preds$lo[1],2),",",round(preds$up[2]-preds$up[1],2),")")
paste0(round(preds$rr,2)," (",round(preds$rr.lo,2),",",round(preds$rr.up,2),")")


######################################
# TABLE 3
######################################

## Regression tables (standardized coefficients)
varz <- c("CONVOYS","COERCION","CONVOYS_PROP","CENT_BTWN_R","DIST_RAIL","URBAN100","OPEN_TERRAIN","OPEN_TERRAIN","ETH_JUD","CONTROL_AREA","DISTLIB_MIN","PEOPLE_PRE","HOUSES_PRE","N_VILLAGES","PEOPLE_KILLED","HOUSES_DESTROYED","REPRISAL")
X <- rayons.panel
X[,which(names(X)%in%varz)] <- scale(X[,which(names(X)%in%varz)])
mod <- gam(PEOPLE_KILLED_p1~as.factor(YRMO)+CONVOYS + COERCION+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN +PEOPLE_PRE+PEOPLE_KILLED+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod1s <- mod
mod <- gam(PEOPLE_KILLED_p1~as.factor(YRMO)+CONVOYS_PROP +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+PEOPLE_PRE+PEOPLE_KILLED+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod2s <- mod
mod <- gam(HOUSES_DESTROYED_p1~as.factor(YRMO)+CONVOYS + COERCION +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+HOUSES_PRE+HOUSES_DESTROYED+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod3s <- mod
mod <- gam(HOUSES_DESTROYED_p1~as.factor(YRMO)+CONVOYS_PROP+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+HOUSES_PRE +HOUSES_DESTROYED+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod4s <- mod
mod <- gam(REPRISAL_p1~as.factor(YRMO)+CONVOYS +COERCION+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+ N_VILLAGES +REPRISAL+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod5s <- mod
mod <- gam(REPRISAL_p1~as.factor(YRMO)+CONVOYS_PROP +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+N_VILLAGES+REPRISAL+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod6s <- mod
varz <- c("CONVOYS","COERCION","CONVOYS_PROP","CENT_BTWN_R","DIST_RAIL","URBAN100","OPEN_TERRAIN","OPEN_TERRAIN","ETH_JUD","CONTROL_AREA","DISTLIB_MIN","PEOPLE_PRE","HOUSES_PRE","N_VILLAGES","PEOPLE_KILLED","HOUSES_DESTROYED","REPRISAL")
covariate.labels <- c("Partisan interdiction","Partisan coercion","Partisan interdiction (proportion)","Rail network centrality","Distance to rail station","Urbanization","Open terrain","Percent Jewish","Partisan control","Distance from frontline","People (pre-war)","People killed (t-1)","Houses (pre-war)","Houses destroyed (t-1)","Villages (pre-war)","Number of reprisals (t-1)")
stargazer(mod1s,mod2s,mod3s,mod4s,mod5s,mod6s,out="Output/TABLE3.tex",star.cutoffs=c(.05,.01,.001),keep = varz,intercept.bottom=T,keep.stat = c("n","ll","adj.rsq","ubre"),covariate.labels=covariate.labels,dep.var.labels =c("People killed","Houses destroyed","Number of reprisals"))

# AIC & EDF
print.xtable(
  xtable(t(data.frame(
    QAIC=r1(c(qaic(mod1s),
              qaic(mod2s),
              qaic(mod3s),
              qaic(mod4s),
              qaic(mod5s),
              qaic(mod6s))),
    DevExp=r3(c(summary(mod1s)$dev.expl,
                summary(mod2s)$dev.expl,
                summary(mod3s)$dev.expl,
                summary(mod4s)$dev.expl,
                summary(mod5s)$dev.expl,
                summary(mod6s)$dev.expl)),
    EDF=r3(c(summary(mod1s)$s.table[1],
             summary(mod2s)$s.table[1],
             summary(mod3s)$s.table[1],
             summary(mod4s)$s.table[1],
             summary(mod5s)$s.table[1],
             summary(mod6s)$s.table[1]))
  ))),
  file="Output/TABLE3_AIC.tex")





######################################
# FIGURE 4
######################################


# Load & crop data
load("Shapes/BELARUS_RAYONS_v8.RData")
borders <- readShapePoly("Shapes/April_30_1941")
proj4string(borders) <- CRS("+proj=aea +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs")
proj4string(map) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
borders <- spTransform(borders,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
borders <- gBuffer(borders, byid=TRUE, width=0)
bb <- bbox(map)
bb[1,1] <- bb[1,1]-(bb[1,2]-bb[1,1])/8
bb[1,2] <- bb[1,2]+(bb[1,2]-bb[1,1])/8
bb[2,1] <- bb[2,1]-(bb[2,2]-bb[2,1])/8
bb[2,2] <- bb[2,2]+(bb[2,2]-bb[2,1])/8
WKTs <- paste("POLYGON((",bb[1,1],bb[2,1],",",bb[1,2],bb[2,1],",",bb[1,2],bb[2,2],",",bb[1,1],bb[2,2],",",bb[1,1],bb[2,1],"))")
test <- readWKT(WKTs)
r <- borders
proj4string(test) <- CRS(proj4string(map))
proj4string(r) <- CRS(proj4string(map))
borders2 <- crop(r, test)

# People killed
p1 <- plot.gam(mod1,n2=100)[[1]]
png("Output/FIG4a.png",width=4,height=3,res=300,units="in")
par(mar=c(0,0,0,0))
plot(test,col="lightblue",border=NA)
plot(borders2,add=T,col="beige",border="grey")
image(p1$x, p1$y, matrix(c(p1$fit),length(p1$y),length(p1$x)), col=rev(gray.colors(256,start=.5)), xlab="",ylab="", axes=FALSE,add=T)
plot(map,add=T,col=NA,border=rgb(.8,.8,.8))
plot(mod1,se=F,shade=T,add=T,lwd=2,n2=100,contour.col="purple")
lls <- c(bbox(map)[1,1]-2.5,mean(c(bbox(map)[2,1]+1.2,1.1*(bbox(map)[2,2]-.7))))
legend(x=bb[1,2]-7.2,y=bb[2,2]+.1,lwd=2,legend="Fitted levels",bty="n",cex=1.2)
dev.off()

# Houses destroyed
p1 <- plot.gam(mod3,n2=100)[[1]]
png("Output/FIG4b.png",width=4,height=3,res=300,units="in")
par(mar=c(0,0,0,0))
plot(test,col="lightblue",border=NA)
plot(borders2,add=T,col="beige",border="grey")
image(p1$x, p1$y, matrix(c(p1$fit),length(p1$y),length(p1$x)), col=rev(gray.colors(256,start=.5)), xlab="",ylab="", axes=FALSE,add=T)
plot(map,add=T,col=NA,border=rgb(.8,.8,.8))
plot(mod3,se=F,shade=T,add=T,lwd=2,n2=100,contour.col="purple")
lls <- c(bbox(map)[1,1]-2.5,mean(c(bbox(map)[2,1]+1.2,1.1*(bbox(map)[2,2]-.7))))
legend(x=bb[1,2]-7.2,y=bb[2,2]+.1,lwd=2,legend="Fitted levels",bty="n",cex=1.2)
dev.off()

# Reprisals
p1 <- plot.gam(mod5,n2=100)[[1]]
png("Output/FIG4c.png",width=4,height=3,res=300,units="in")
par(mar=c(0,0,0,0))
plot(test,col="lightblue",border=NA)
plot(borders2,add=T,col="beige",border="grey")
image(p1$x, p1$y, matrix(c(p1$fit),length(p1$y),length(p1$x)), col=rev(gray.colors(256,start=.5)), xlab="",ylab="", axes=FALSE,add=T)
plot(map,add=T,col=NA,border=rgb(.8,.8,.8))
plot(mod5,se=F,shade=T,add=T,lwd=2,n2=100,contour.col="purple")
lls <- c(bbox(map)[1,1]-2.5,mean(c(bbox(map)[2,1]+1.2,1.1*(bbox(map)[2,2]-.7))))
legend(x=bb[1,2]-7.2,y=bb[2,2]+.1,lwd=2,legend="Fitted levels",bty="n",cex=1.2)
dev.off()


######################################
# FIGURE 5
######################################

# People killed
disag <- as.character(sort(unique(mod1$model$`as.factor(YRMO)`)))
preds.t <- lapply(1:length(disag),function(t){preds <- predmanGAM(mod=mod1,var="CONVOYS",vals=c(0),quadratic=F,fax=list(YRMO=disag[t]))
data.frame(y=preds$y,lo=preds$lo,up=preds$up)
})
preds.t <- do.call(rbind,preds.t)
preds.t$YRMO <- disag; preds.t$TID <- 1:length(disag)
preds.t1 <- preds.t
# Houses destroyed
disag <- as.character(sort(unique(mod3$model$`as.factor(YRMO)`)))
preds.t <- lapply(1:length(disag),function(t){preds <- predmanGAM(mod=mod3,var="CONVOYS",vals=c(0),quadratic=F,fax=list(YRMO=disag[t]))
data.frame(y=preds$y,lo=preds$lo,up=preds$up)
})
preds.t <- do.call(rbind,preds.t)
preds.t$YRMO <- disag; preds.t$TID <- 1:length(disag)
preds.t2 <- preds.t
# Number of reprisals
disag <- as.character(sort(unique(mod5$model$`as.factor(YRMO)`)))
preds.t <- lapply(1:length(disag),function(t){preds <- predmanGAM(mod=mod5,var="CONVOYS",vals=c(0),quadratic=F,fax=list(YRMO=disag[t]))
data.frame(y=preds$y,lo=preds$lo,up=preds$up)
})
preds.t <- do.call(rbind,preds.t)
preds.t$YRMO <- disag; preds.t$TID <- 1:length(disag)
preds.t3 <- preds.t


# plot
ymaxx <- 200
ymaxx2 <- 3
png("Output/FIG5.png",width=8,height=2,res=300,units="in")
par(mar=c(2,2,0.5,2))
plot(0,0,col=NA,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",xlim=range(preds.t$TID),ylim=c(0,ymaxx))
axis(1,at=seq(2,max(preds.t$TID),by=6),labels=paste0(substr(preds.t$YRMO,5,6),"/",substr(preds.t$YRMO,1,4))[seq(2,max(preds.t$TID),by=6)],cex.axis=.75,lwd=2)
axis(1,at=1:max(preds.t$TID),labels=FALSE,cex.axis=.5)
axis(2,at=seq(0,ymaxx,by=ymaxx/5),labels=seq(0,ymaxx,by=ymaxx/5),cex.axis=.75,lwd=2)
axis(2,at=seq(0,ymaxx,by=ymaxx/20),labels=FALSE)
segments(x0=preds.t1$TID-.3,x1=preds.t1$TID-.3,y0=preds.t1$lo,y1=preds.t1$up,col=rgb(.9,.6,.6),lwd=3,lend=4)
points(x=preds.t1$TID-.3,y=preds.t1$y,pch=16,col="darkred",cex=.5)
segments(x0=preds.t2$TID,x1=preds.t2$TID,y0=preds.t2$lo,y1=preds.t2$up,col=rgb(.6,.6,.9),lwd=3,lend=4)
points(x=preds.t2$TID,y=preds.t2$y,pch=16,col="darkblue",cex=.5)
par(new=T)
plot(0,0,col=NA,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",xlim=range(preds.t$TID),ylim=c(0,ymaxx2))
axis(4,at=seq(0,ymaxx2,by=ymaxx2/3),labels=seq(0,ymaxx2,by=ymaxx2/3),cex.axis=.75,lwd=2,col="darkgreen",col.axis="darkgreen")
axis(4,at=seq(0,ymaxx2,by=ymaxx2/15),labels=FALSE,col="darkgreen",col.axis="darkgreen")
segments(x0=preds.t3$TID+.3,x1=preds.t3$TID+.3,y0=preds.t3$lo,y1=preds.t3$up,col=rgb(.6,.9,.6),lwd=3,lend=4)
points(x=preds.t3$TID+.3,y=preds.t3$y,pch=16,col="darkgreen",cex=.5)
legend("topleft",fill=c("darkred","darkblue"),legend=c("Houses destroyed","People killed"),bty="n")
legend("topright",fill=c("darkgreen"),legend=c("Number of\nreprisals"),bty="n")
dev.off()


# plot (b/w)
ymaxx <- 200
ymaxx2 <- 3
#png("Output/FIG5_bw.png",width=8,height=2,res=300,units="in")
tiff("Output/FIG5_bw.tiff",width=8,height=2,res=300,units="in")
layout(matrix(c(rep(1,7),2),ncol=1))
par(mar=c(2,2,0.5,2))
plot(0,0,col=NA,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",xlim=range(preds.t$TID),ylim=c(0,ymaxx))
axis(1,at=seq(2,max(preds.t$TID),by=6),labels=paste0(substr(preds.t$YRMO,5,6),"/",substr(preds.t$YRMO,1,4))[seq(2,max(preds.t$TID),by=6)],cex.axis=.75,lwd=2)
axis(1,at=1:max(preds.t$TID),labels=FALSE,cex.axis=.5)
axis(2,at=seq(0,ymaxx,by=ymaxx/5),labels=seq(0,ymaxx,by=ymaxx/5),cex.axis=.75,lwd=2)
axis(2,at=seq(0,ymaxx,by=ymaxx/20),labels=FALSE)
segments(x0=preds.t1$TID-.3,x1=preds.t1$TID-.3,y0=preds.t1$lo,y1=preds.t1$up,col="gray20",lwd=2,lend=4,lty=1)
points(x=preds.t1$TID-.3,y=preds.t1$y,pch=1,col="gray50",cex=.5)
segments(x0=preds.t2$TID,x1=preds.t2$TID,y0=preds.t2$lo,y1=preds.t2$up,col="gray50",lwd=2,lend=4,lty=1)
points(x=preds.t2$TID,y=preds.t2$y,pch=2,col="black",cex=.5)
par(new=T)
plot(0,0,col=NA,bty="n",xaxt="n",yaxt="n",xlab="",ylab="",xlim=range(preds.t$TID),ylim=c(0,ymaxx2))
axis(4,at=seq(0,ymaxx2,by=ymaxx2/3),labels=seq(0,ymaxx2,by=ymaxx2/3),cex.axis=.75,lwd=2,col="black",col.axis="black")
axis(4,at=seq(0,ymaxx2,by=ymaxx2/15),labels=FALSE,col="black",col.axis="black")
segments(x0=preds.t3$TID+.3,x1=preds.t3$TID+.3,y0=preds.t3$lo,y1=preds.t3$up,col="gray80",lwd=2,lend=4)
points(x=preds.t3$TID+.3,y=preds.t3$y,pch=3,col="black",cex=.5)
par(mar=c(0,0,0,0))
plot.new()
legend("top",legend=c("Houses destroyed","People killed","Number of reprisals (right axis)"),bty="n",pch=c(1:3),ncol=3,col="gray40",text.col = "gray40")
dev.off()



## Distance to Germany
head(map)
temp <- data.frame(RAYON_ID=map$RAYON_ID,DIST2BERLIN=NA)
berlin <- c(13.383333,52.516667)
coords <- coordinates(map)
for(i in 1:nrow(coords)){temp[i,"DIST2BERLIN"] <- dist(rbind(berlin,coords[i,]))}


## Regression tables (standardized coefficients)
varz <- c("CONVOYS","COERCION","CONVOYS_PROP","CENT_BTWN_R","DIST_RAIL","URBAN100","OPEN_TERRAIN","OPEN_TERRAIN","ETH_JUD","CONTROL_AREA","DISTLIB_MIN","PEOPLE_PRE","HOUSES_PRE","N_VILLAGES","PEOPLE_KILLED","HOUSES_DESTROYED","REPRISAL","DIST2BERLIN")
X <- rayons.panel
X <- merge(X,temp,by="RAYON_ID",all.x=T,all.y=F)
X[,which(names(X)%in%varz)] <- scale(X[,which(names(X)%in%varz)])
mod <- gam(PEOPLE_KILLED_p1~as.factor(YRMO)+CONVOYS + COERCION+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN +PEOPLE_PRE+PEOPLE_KILLED+DIST2BERLIN+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod1s <- mod
mod <- gam(PEOPLE_KILLED_p1~as.factor(YRMO)+CONVOYS_PROP +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+PEOPLE_PRE+PEOPLE_KILLED+DIST2BERLIN+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod2s <- mod
mod <- gam(HOUSES_DESTROYED_p1~as.factor(YRMO)+CONVOYS + COERCION +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+HOUSES_PRE+HOUSES_DESTROYED+DIST2BERLIN+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod3s <- mod
mod <- gam(HOUSES_DESTROYED_p1~as.factor(YRMO)+CONVOYS_PROP+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+HOUSES_PRE +HOUSES_DESTROYED+DIST2BERLIN+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod4s <- mod
mod <- gam(REPRISAL_p1~as.factor(YRMO)+CONVOYS +COERCION+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+ N_VILLAGES +REPRISAL+DIST2BERLIN+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod5s <- mod
mod <- gam(REPRISAL_p1~as.factor(YRMO)+CONVOYS_PROP +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+N_VILLAGES+REPRISAL+DIST2BERLIN+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod6s <- mod

varz <- c("CONVOYS","COERCION","CONVOYS_PROP","CENT_BTWN_R","DIST_RAIL","URBAN100","OPEN_TERRAIN","OPEN_TERRAIN","ETH_JUD","CONTROL_AREA","DISTLIB_MIN","PEOPLE_PRE","HOUSES_PRE","N_VILLAGES","PEOPLE_KILLED","HOUSES_DESTROYED","REPRISAL","DIST2BERLIN")
covariate.labels <- c("Partisan interdiction","Partisan coercion","Partisan interdiction (proportion)","Rail network centrality","Distance to rail station","Urbanization","Open terrain","Percent Jewish","Partisan control","Distance from frontline","People (pre-war)","People killed (t-1)","Houses (pre-war)","Houses destroyed (t-1)","Villages (pre-war)","Number of reprisals (t-1)","Distance to Berlin")
stargazer(mod1s,mod2s,mod3s,mod4s,mod5s,mod6s,out="Output/TABLE3_berlin.tex",star.cutoffs=c(.05,.01,.001),keep = varz,intercept.bottom=T,keep.stat = c("n","ll","adj.rsq","ubre"),covariate.labels=covariate.labels,dep.var.labels =c("People killed","Houses destroyed","Number of reprisals"))

# AIC & EDF
print.xtable(
  xtable(t(data.frame(
    QAIC=r1(c(qaic(mod1s),
              qaic(mod2s),
              qaic(mod3s),
              qaic(mod4s),
              qaic(mod5s),
              qaic(mod6s))),
    DevExp=r3(c(summary(mod1s)$dev.expl,
                summary(mod2s)$dev.expl,
                summary(mod3s)$dev.expl,
                summary(mod4s)$dev.expl,
                summary(mod5s)$dev.expl,
                summary(mod6s)$dev.expl)),
    EDF=r3(c(summary(mod1s)$s.table[1],
             summary(mod2s)$s.table[1],
             summary(mod3s)$s.table[1],
             summary(mod4s)$s.table[1],
             summary(mod5s)$s.table[1],
             summary(mod6s)$s.table[1]))
  ))),
  file="Output/TABLE3_berlin_AIC.tex")




## Distance to Germany (no spline)
head(map)
temp <- data.frame(RAYON_ID=map$RAYON_ID,DIST2BERLIN=NA)
berlin <- c(13.383333,52.516667)
coords <- coordinates(map)
for(i in 1:nrow(coords)){temp[i,"DIST2BERLIN"] <- dist(rbind(berlin,coords[i,]))}


## Regression tables (standardized coefficients)
varz <- c("CONVOYS","COERCION","CONVOYS_PROP","CENT_BTWN_R","DIST_RAIL","URBAN100","OPEN_TERRAIN","OPEN_TERRAIN","ETH_JUD","CONTROL_AREA","DISTLIB_MIN","PEOPLE_PRE","HOUSES_PRE","N_VILLAGES","PEOPLE_KILLED","HOUSES_DESTROYED","REPRISAL","DIST2BERLIN")
X <- rayons.panel
X <- merge(X,temp,by="RAYON_ID",all.x=T,all.y=F)
X[,which(names(X)%in%varz)] <- scale(X[,which(names(X)%in%varz)])
mod <- gam(PEOPLE_KILLED_p1~as.factor(YRMO)+CONVOYS + COERCION+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN +PEOPLE_PRE+PEOPLE_KILLED+DIST2BERLIN,data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod1s <- mod
mod <- gam(PEOPLE_KILLED_p1~as.factor(YRMO)+CONVOYS_PROP +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+PEOPLE_PRE+PEOPLE_KILLED+DIST2BERLIN,data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod2s <- mod
mod <- gam(HOUSES_DESTROYED_p1~as.factor(YRMO)+CONVOYS + COERCION +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+HOUSES_PRE+HOUSES_DESTROYED+DIST2BERLIN,data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod3s <- mod
mod <- gam(HOUSES_DESTROYED_p1~as.factor(YRMO)+CONVOYS_PROP+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+HOUSES_PRE +HOUSES_DESTROYED+DIST2BERLIN,data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod4s <- mod
mod <- gam(REPRISAL_p1~as.factor(YRMO)+CONVOYS +COERCION+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+ N_VILLAGES +REPRISAL+DIST2BERLIN,data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod5s <- mod
mod <- gam(REPRISAL_p1~as.factor(YRMO)+CONVOYS_PROP +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+N_VILLAGES+REPRISAL+DIST2BERLIN,data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod6s <- mod

varz <- c("CONVOYS","COERCION","CONVOYS_PROP","CENT_BTWN_R","DIST_RAIL","URBAN100","OPEN_TERRAIN","OPEN_TERRAIN","ETH_JUD","CONTROL_AREA","DISTLIB_MIN","PEOPLE_PRE","HOUSES_PRE","N_VILLAGES","PEOPLE_KILLED","HOUSES_DESTROYED","REPRISAL","DIST2BERLIN")
covariate.labels <- c("Partisan interdiction","Partisan coercion","Partisan interdiction (proportion)","Rail network centrality","Distance to rail station","Urbanization","Open terrain","Percent Jewish","Partisan control","Distance from frontline","People (pre-war)","People killed (t-1)","Houses (pre-war)","Houses destroyed (t-1)","Villages (pre-war)","Number of reprisals (t-1)","Distance to Berlin")
stargazer(mod1s,mod2s,mod3s,mod4s,mod5s,mod6s,out="Output/TABLE3_berlinns.tex",star.cutoffs=c(.05,.01,.001),keep = varz,intercept.bottom=T,keep.stat = c("n","ll","adj.rsq","ubre"),covariate.labels=covariate.labels,dep.var.labels =c("People killed","Houses destroyed","Number of reprisals"))

# AIC & EDF
print.xtable(
  xtable(t(data.frame(
    QAIC=r1(c(qaic(mod1s),
              qaic(mod2s),
              qaic(mod3s),
              qaic(mod4s),
              qaic(mod5s),
              qaic(mod6s))),
    DevExp=r3(c(summary(mod1s)$dev.expl,
                summary(mod2s)$dev.expl,
                summary(mod3s)$dev.expl,
                summary(mod4s)$dev.expl,
                summary(mod5s)$dev.expl,
                summary(mod6s)$dev.expl))
  )
  )),
  file="Output/TABLE3_berlinns_AIC.tex")







######################################
# TABLE 4
######################################

## Regression tables (standardized coefficients)
varz <- c("CONVOYS","COERCION","CONVOYS_PROP","CENT_BTWN_R","DIST_RAIL","URBAN100","OPEN_TERRAIN","OPEN_TERRAIN","ETH_JUD","CONTROL_AREA","DISTLIB_MIN","PEOPLE_PRE","HOUSES_PRE","N_VILLAGES","PEOPLE_KILLED","HOUSES_DESTROYED","REPRISAL")
X <- villages.panel
X[,which(names(X)%in%varz)] <- scale(X[,which(names(X)%in%varz)])
mod <- gam(PEOPLE_KILLED_p1~as.factor(YRMO)+CONVOYS + COERCION+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN +PEOPLE_PRE+PEOPLE_KILLED+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod1s <- mod
mod <- gam(PEOPLE_KILLED_p1~as.factor(YRMO)+CONVOYS_PROP +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+PEOPLE_PRE+PEOPLE_KILLED+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod2s <- mod
mod <- gam(HOUSES_DESTROYED_p1~as.factor(YRMO)+CONVOYS + COERCION +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+HOUSES_PRE+HOUSES_DESTROYED+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod3s <- mod
mod <- gam(HOUSES_DESTROYED_p1~as.factor(YRMO)+CONVOYS_PROP+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+HOUSES_PRE +HOUSES_DESTROYED+s(LONG,LAT),data=X,family="quasipoisson");summary(mod);#qaic(mod)
mod4s <- mod
mod <- gam(REPRISAL_p1~as.factor(YRMO)+CONVOYS +COERCION+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+ PEOPLE_PRE +REPRISAL+s(LONG,LAT),data=X,family="logit");summary(mod);#qaic(mod)
mod5s <- mod
mod <- gam(REPRISAL_p1~as.factor(YRMO)+CONVOYS_PROP +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+PEOPLE_PRE+REPRISAL+s(LONG,LAT),data=X,family="logit");summary(mod);#qaic(mod)
mod6s <- mod

varz <- c("CONVOYS","COERCION","CONVOYS_PROP","CENT_BTWN_R","DIST_RAIL","URBAN100","OPEN_TERRAIN","OPEN_TERRAIN","ETH_JUD","CONTROL_AREA","DISTLIB_MIN","PEOPLE_PRE","HOUSES_PRE","PEOPLE_KILLED","HOUSES_DESTROYED","REPRISAL")
covariate.labels <- c("Partisan interdiction","Partisan coercion","Partisan interdiction (proportion)","Rail network centrality","Distance to rail station","Urbanization","Open terrain","Percent Jewish","Partisan control","Distance from frontline","People (pre-war)","People killed (t-1)","Houses (pre-war)","Houses destroyed (t-1)","Number of reprisals (t-1)")
stargazer(mod1s,mod2s,mod3s,mod4s,mod5s,mod6s,out="Output/TABLE5.tex",star.cutoffs=c(.05,.01,.001),keep = varz,intercept.bottom=T,keep.stat = c("n","ll","adj.rsq","ubre"),covariate.labels=covariate.labels,dep.var.labels =c("People killed","Houses destroyed","Number of reprisals"))

# AIC & EDF
print.xtable(
  xtable(t(data.frame(
    #     QAIC=r1(c(qaic(mod1s),
    #               qaic(mod2s),
    #               qaic(mod3s),
    #               qaic(mod4s),
    #               qaic(mod5s),
    #               qaic(mod6s))),
    DevExp=r3(c(summary(mod1s)$dev.expl,
                summary(mod2s)$dev.expl,
                summary(mod3s)$dev.expl,
                summary(mod4s)$dev.expl,
                summary(mod5s)$dev.expl,
                summary(mod6s)$dev.expl)),
    EDF=r3(c(summary(mod1s)$s.table[1],
             summary(mod2s)$s.table[1],
             summary(mod3s)$s.table[1],
             summary(mod4s)$s.table[1],
             summary(mod5s)$s.table[1],
             summary(mod6s)$s.table[1]))
  ))),
  file="Output/TABLE5_AIC.tex")



######################################
# FIGURE 3b
######################################

load("Data/BELARUS_VILLAGES_MONTHS_Full.RData")
head(villages.panel)

mod <- gam(HOUSES_DESTROYED_p1~as.factor(YRMO)+CONVOYS + COERCION +CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+HOUSES_PRE+HOUSES_DESTROYED+s(LONG,LAT),data=villages.panel,family="quasipoisson");summary(mod);#qaic(mod)
mod3 <- mod; rm(mod)
#plot(mod3)

mod <- gam(HOUSES_DESTROYED_p1~as.factor(YRMO)+CONVOYS_PROP+CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN+HOUSES_PRE +HOUSES_DESTROYED+s(LONG,LAT),data=villages.panel,family="quasipoisson");summary(mod);#qaic(mod)
mod4 <- mod; rm(mod)

# CF plots
png("Output/FIG3bs.png",height=3,width=12,units="in",res=300)
par(mfrow=c(1,3),mar=c(4,4,.5,.5))
mod <- mod3
contplot_GAM(mod,var="CONVOYS",xlab="Partisan interdiction (convoys derailed)",ylab="Houses destroyed by Germans",cax=1.4,fax=list(YRMO=194310, VID=118),col.ci=rgb(.5,.33,1,.75))
contplot_GAM(mod,var="COERCION",xlab="Partisan coercion (garrisons attacked)",ylab="Houses destroyed by Germans",cax=1.4,fax=list(YRMO=194310, VID=118))
mod <- mod4
contplot_GAM(mod,var="CONVOYS_PROP",xlab="Partisan interdiction (proportion)",ylab="Houses destroyed by Germans",cax=1.4,fax=list(YRMO=194310, VID=118),col.ci=rgb(.5,.33,1,.75))
dev.off()

rm(mod,mod3,mod4)


## CFs
# People killed
preds <- predmanGAM(mod=mod1,var="CONVOYS",vals=c(0,10),quadratic=F,fax=list(YRMO=194310))
paste0(round(preds$y[2]-preds$y[1],2)," (",round(preds$lo[2]-preds$lo[1],2),",",round(preds$up[2]-preds$up[1],2),")")
paste0(round(preds$rr,2)," (",round(preds$rr.lo,2),",",round(preds$rr.up,2),")")

preds <- predmanGAM(mod=mod1,var="COERCION",vals=c(0,1),quadratic=F,fax=list(YRMO=194310))
paste0(round(preds$y[2]-preds$y[1],2)," (",round(preds$lo[2]-preds$lo[1],2),",",round(preds$up[2]-preds$up[1],2),")")
paste0(round(preds$rr,2)," (",round(preds$rr.lo,2),",",round(preds$rr.up,2),")")

preds <- predmanGAM(mod=mod2,var="CONVOYS_PROP",vals=c(0,1),quadratic=F,fax=list(YRMO=194310))
paste0(round(preds$y[2]-preds$y[1],2)," (",round(preds$lo[2]-preds$lo[1],2),",",round(preds$up[2]-preds$up[1],2),")")
paste0(round(preds$rr,2)," (",round(preds$rr.lo,2),",",round(preds$rr.up,2),")")

preds <- predmanGAM(mod=mod1,var="CONTROL_AREA",vals=c(0,.5,1),quadratic=F,fax=list(YRMO=194310))
paste0(round(preds$y[1],2)," (",round(preds$lo[1],2),",",round(preds$up[1],2),")")
paste0(round(preds$y[3],2)," (",round(preds$lo[3],2),",",round(preds$up[3],2),")")
paste0(round(preds$y[2],2)," (",round(preds$lo[2],2),",",round(preds$up[2],2),")")
paste0(round(preds$y[2]-preds$y[1],2)," (",round(preds$lo[2]-preds$lo[1],2),",",round(preds$up[2]-preds$up[1],2),")")
paste0(round(preds$rr,2)," (",round(preds$rr.lo,2),",",round(preds$rr.up,2),")")

preds <- predmanGAM(mod=mod1,var="ETH_JUD",vals=c(1,30),quadratic=F,fax=list(YRMO=194310))
paste0(round(preds$y[2]-preds$y[1],2)," (",round(preds$lo[2]-preds$lo[1],2),",",round(preds$up[2]-preds$up[1],2),")")
paste0(round(preds$rr,2)," (",round(preds$rr.lo,2),",",round(preds$rr.up,2),")")

summary(villages.panel$DISTLIB_MIN)
preds <- predmanGAM(mod=mod1,var="DISTLIB_MIN",vals=c(77,578),quadratic=F,fax=list(YRMO=194310))
paste0(round(preds$y[2]-preds$y[1],2)," (",round(preds$lo[2]-preds$lo[1],2),",",round(preds$up[2]-preds$up[1],2),")")
paste0(round(preds$rr,2)," (",round(preds$rr.lo,2),",",round(preds$rr.up,2),")")




######################################
# TABLE 5
######################################

## Rayon-level
load("Data/BELARUS_RAYONS_MONTHS_Full.RData")
load("Data/BELARUS_RAYONS_Full_sp.RData")
map_crd <- data.frame(RAYON_ID=map$RAYON_ID,LONG=coordinates(map)[,1],LAT=coordinates(map)[,2])
rayons.panel <- merge(rayons.panel,map_crd,by="RAYON_ID")
rayons.panel <- merge(rayons.panel,data.frame(YRMO=sort(unique(rayons.panel$YRMO)),TID=1:length(unique(rayons.panel$YRMO))),by="YRMO")
X <- rayons.panel
drange <- 194308:194310
disag <- unique(map$RAYON_ID)

# Pre/post variables
i <- 23
dim(X)
mat.list <- lapply(1:length(disag),function(i){#print(paste(tvar,i))
  sub. <- X[X$RAYON_ID==disag[i],]
  tidz <- unique(sub.$TID[sub.$YRMO%in%drange])
  sub. <- sub.[order(sub.$TID),]
  sub. <- sub.[!duplicated(sub.$TID,fromLast=T),]
  sub2 <- data.frame(RAYON_ID=disag[i])
  sub2$P_PRE <- sum(sub.$PEOPLE_KILLED[sub.$TID<min(tidz)])
  sub2$P_PRE1 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$P_PRE2 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$P_PRE3 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$P_PRE4 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$P_PRE5 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$P_PRE6 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$P_DURING <- sum(sub.$PEOPLE_KILLED[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$P_POST <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)])
  sub2$P_POST1 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+1)])
  sub2$P_POST2 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+2)])
  sub2$P_POST3 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+3)])
  sub2$P_POST4 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+4)])
  sub2$P_POST5 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+5)])
  sub2$P_POST6 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+6)])
  sub2$H_PRE <- sum(sub.$HOUSES_DESTROYED[sub.$TID<min(tidz)])
  sub2$H_PRE1 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$H_PRE2 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$H_PRE3 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$H_PRE4 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$H_PRE5 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$H_PRE6 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$H_DURING <- sum(sub.$HOUSES_DESTROYED[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$H_POST <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)])
  sub2$H_POST1 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+1)])
  sub2$H_POST2 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+2)])
  sub2$H_POST3 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+3)])
  sub2$H_POST4 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+4)])
  sub2$H_POST5 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+5)])
  sub2$H_POST6 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+6)])
  sub2$R_PRE <- sum(sub.$REPRISAL[sub.$TID<min(tidz)])
  sub2$R_PRE1 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$R_PRE2 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$R_PRE3 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$R_PRE4 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$R_PRE5 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$R_PRE6 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$R_DURING <- sum(sub.$REPRISAL[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$R_POST <- sum(sub.$REPRISAL[sub.$TID>max(tidz)])
  sub2$R_POST1 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+1)])
  sub2$R_POST2 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+2)])
  sub2$R_POST3 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+3)])
  sub2$R_POST4 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+4)])
  sub2$R_POST5 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+5)])
  sub2$R_POST6 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+6)])
  sub2$T_CONVOYS_PRE <- sum(sub.$CONVOYS[sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE1 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE2 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE3 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE4 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE5 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE6 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_DURING <- sum(sub.$CONVOYS[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$T_COERCION_PRE <- sum(sub.$COERCION[sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE1 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE2 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE3 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE4 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE5 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE6 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$T_COERCION_DURING <- sum(sub.$COERCION[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$DISTLIB_MIN_PRE <- mean(sub.$DISTLIB_MIN[sub.$TID==min(tidz)])
  sub2
})
mat <- do.call(rbind,mat.list)
head(mat)

# Diff-in-diff
mat$T_CONVOYS <- 1*(mat$T_CONVOYS_DURING>mean(mat$T_CONVOYS_DURING))
mat$T_COERCION <- 1*(mat$T_COERCION_DURING>mean(mat$T_COERCION_DURING))
mat$ACTIVE <- 1*(mat$T_COERCION+mat$T_CONVOYS>0)
mat_pre <- mat
mat_pre$POST <- 0
mat_pre$P_ALL <- mat_pre$P_PRE
mat_pre$P_ALL3 <- mat_pre$P_PRE3
mat_pre$H_ALL <- mat_pre$H_PRE
mat_pre$H_ALL3 <- mat_pre$H_PRE3
mat_pre$R_ALL <- mat_pre$R_PRE
mat_pre$R_ALL3 <- mat_pre$R_PRE3
mat_pre$P_DALL <- mat_pre$P_PRE
mat_pre$P_DALL3 <- mat_pre$P_PRE3
mat_pre$H_DALL <- mat_pre$H_PRE
mat_pre$H_DALL3 <- mat_pre$H_PRE3
mat_pre$R_DALL <- mat_pre$R_PRE
mat_pre$R_DALL3 <- mat_pre$R_PRE3
mat_post <- mat
mat_post$POST <- 1
mat_post$P_ALL <- mat_post$P_POST
mat_post$P_ALL3 <- mat_post$P_POST3
mat_post$H_ALL <- mat_post$H_POST
mat_post$H_ALL3 <- mat_post$H_POST3
mat_post$R_ALL <- mat_post$R_POST
mat_post$R_ALL3 <- mat_post$R_POST3
mat_post$P_DALL <- mat_post$P_POST+mat_post$P_DURING
mat_post$P_DALL3 <- mat_post$P_POST3+mat_post$P_DURING
mat_post$H_DALL <- mat_post$H_POST+mat_post$P_DURING
mat_post$H_DALL3 <- mat_post$H_POST3+mat_post$P_DURING
mat_post$R_DALL <- mat_post$R_POST+mat_post$P_DURING
mat_post$R_DALL3 <- mat_post$R_POST3+mat_post$P_DURING
mat_dind <- rbind(mat_pre,mat_post)
names(mat_dind)


## Village-level
load("Data/BELARUS_VILLAGES_MONTHS_Full.RData")
load("Data/BELARUS_VILLAGES_Full.RData")
villages.panel <- merge(villages.panel,data.frame(YRMO=sort(unique(villages.panel$YRMO)),TID=1:length(unique(villages.panel$YRMO))),by="YRMO")
X <- villages.panel
drange <- 194308:194310
disag <- unique(villages$VID)

# Pre/post variables
i <- 23
dim(X)
mat.list <- lapply(1:length(disag),function(i){print(i)
  sub. <- X[X$VID==disag[i],]
  tidz <- unique(sub.$TID[sub.$YRMO%in%drange])
  sub. <- sub.[order(sub.$TID),]
  sub. <- sub.[!duplicated(sub.$TID,fromLast=T),]
  sub2 <- data.frame(VID=disag[i])
  sub2$P_PRE <- sum(sub.$PEOPLE_KILLED[sub.$TID<min(tidz)])
  sub2$P_PRE1 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$P_PRE2 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$P_PRE3 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$P_PRE4 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$P_PRE5 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$P_PRE6 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$P_DURING <- sum(sub.$PEOPLE_KILLED[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$P_POST <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)])
  sub2$P_POST1 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+1)])
  sub2$P_POST2 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+2)])
  sub2$P_POST3 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+3)])
  sub2$P_POST4 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+4)])
  sub2$P_POST5 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+5)])
  sub2$P_POST6 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+6)])
  sub2$H_PRE <- sum(sub.$HOUSES_DESTROYED[sub.$TID<min(tidz)])
  sub2$H_PRE1 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$H_PRE2 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$H_PRE3 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$H_PRE4 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$H_PRE5 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$H_PRE6 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$H_DURING <- sum(sub.$HOUSES_DESTROYED[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$H_POST <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)])
  sub2$H_POST1 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+1)])
  sub2$H_POST2 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+2)])
  sub2$H_POST3 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+3)])
  sub2$H_POST4 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+4)])
  sub2$H_POST5 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+5)])
  sub2$H_POST6 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+6)])
  sub2$R_PRE <- sum(sub.$REPRISAL[sub.$TID<min(tidz)])
  sub2$R_PRE1 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$R_PRE2 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$R_PRE3 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$R_PRE4 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$R_PRE5 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$R_PRE6 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$R_DURING <- sum(sub.$REPRISAL[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$R_POST <- sum(sub.$REPRISAL[sub.$TID>max(tidz)])
  sub2$R_POST1 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+1)])
  sub2$R_POST2 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+2)])
  sub2$R_POST3 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+3)])
  sub2$R_POST4 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+4)])
  sub2$R_POST5 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+5)])
  sub2$R_POST6 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+6)])
  sub2$T_CONVOYS_PRE <- sum(sub.$CONVOYS[sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE1 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE2 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE3 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE4 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE5 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE6 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_DURING <- sum(sub.$CONVOYS[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$T_COERCION_PRE <- sum(sub.$COERCION[sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE1 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE2 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE3 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE4 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE5 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE6 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$T_COERCION_DURING <- sum(sub.$COERCION[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$DISTLIB_MIN_PRE <- mean(sub.$DISTLIB_MIN[sub.$TID==min(tidz)])
  sub2
})
mat <- do.call(rbind,mat.list)
head(mat)
save(mat, file="Data/BELARUS_VILLAGES_dindmat.RData")
load("Data/BELARUS_VILLAGES_dindmat.RData")
sum(mat$T_CONVOYS_DURING>0)

# Diff-in-diff
mat$T_CONVOYS <- 1*(mat$T_CONVOYS_DURING>mean(mat$T_CONVOYS_DURING))
mat$T_COERCION <- 1*(mat$T_COERCION_DURING>mean(mat$T_COERCION_DURING))
mat$T_CONVOYS <- 1*(mat$T_CONVOYS_DURING>0)
mat$T_COERCION <- 1*(mat$T_COERCION_DURING>0)
mat$ACTIVE <- 1*(mat$T_COERCION+mat$T_CONVOYS>0)
mat_pre <- mat
mat_pre$POST <- 0
mat_pre$P_ALL <- mat_pre$P_PRE
mat_pre$P_ALL3 <- mat_pre$P_PRE3
mat_pre$H_ALL <- mat_pre$H_PRE
mat_pre$H_ALL3 <- mat_pre$H_PRE3
mat_pre$R_ALL <- mat_pre$R_PRE
mat_pre$R_ALL3 <- mat_pre$R_PRE3
mat_pre$P_DALL <- mat_pre$P_PRE
mat_pre$P_DALL3 <- mat_pre$P_PRE3
mat_pre$H_DALL <- mat_pre$H_PRE
mat_pre$H_DALL3 <- mat_pre$H_PRE3
mat_pre$R_DALL <- mat_pre$R_PRE
mat_pre$R_DALL3 <- mat_pre$R_PRE3
mat_post <- mat
mat_post$POST <- 1
mat_post$P_ALL <- mat_post$P_POST
mat_post$P_ALL3 <- mat_post$P_POST3
mat_post$H_ALL <- mat_post$H_POST
mat_post$H_ALL3 <- mat_post$H_POST3
mat_post$R_ALL <- mat_post$R_POST
mat_post$R_ALL3 <- mat_post$R_POST3
mat_post$P_DALL <- mat_post$P_POST+mat_post$P_DURING
mat_post$P_DALL3 <- mat_post$P_POST3+mat_post$P_DURING
mat_post$H_DALL <- mat_post$H_POST+mat_post$P_DURING
mat_post$H_DALL3 <- mat_post$H_POST3+mat_post$P_DURING
mat_post$R_DALL <- mat_post$R_POST+mat_post$P_DURING
mat_post$R_DALL3 <- mat_post$R_POST3+mat_post$P_DURING
mat_dind <- rbind(mat_pre,mat_post)
mat_dind2 <- merge(mat_dind,villages,by="VID",all.x=T,all.y=F)
save(mat_dind2, file="Data/BELARUS_VILLAGES_dind.RData")
write.dta(mat_dind2, file="Data/BELARUS_VILLAGES_dind.dta",convert.factors = c( "string"))


## Matching
load("Data/BELARUS_VILLAGES_dind.RData")
X <- mat_dind2[,c("VID","RAYON_ID","POST","T_CONVOYS","T_COERCION","P_ALL","H_ALL","R_ALL","P_ALL3","H_ALL3","R_ALL3","P_DALL","H_DALL","R_DALL","P_DALL3","H_DALL3","R_DALL3","P_DURING","H_DURING","R_DURING","P_PRE","H_PRE","R_PRE","T_CONVOYS_PRE","T_COERCION_PRE","P_PRE3","H_PRE3","R_PRE3","T_CONVOYS_PRE3","T_COERCION_PRE3","CENT_BTWN_R","DIST_RAIL","URBAN100","OPEN_TERRAIN","ETH_JUD","CONTROL_AREA","DISTLIB_MIN_PRE")]

# COERCION
m.out <- matchit(T_COERCION ~ P_PRE + T_COERCION_PRE+ T_CONVOYS_PRE +P_PRE3 + T_COERCION_PRE3 + T_CONVOYS_PRE3  + CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN_PRE, distance="logit",exact=c("RAYON_ID"), ratio=1, data=X, discard="both",reestimate=T,caliper=.25)
summary(m.out,standardize=T)
#plot(m.out)
m.data <- match.data(m.out)
head(m.data)
dim(m.data)
write.dta(m.data, file="Data/BELARUS_VILLAGES_matched_v2.dta",convert.factors = c( "string"))

# CONVOYS
m.out <- matchit(T_CONVOYS ~ P_PRE + T_COERCION_PRE+ T_CONVOYS_PRE +P_PRE3 + T_COERCION_PRE3 + T_CONVOYS_PRE3  + CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN_PRE, distance="logit",exact=c("RAYON_ID"), ratio=1, data=X, discard="both",reestimate=T,caliper=.25)
summary(m.out,standardize=T)
xtable(m.out)
m.data <- match.data(m.out)
head(m.data)
dim(m.data)
write.dta(m.data, file="Data/BELARUS_VILLAGES_matched_v1.dta",convert.factors = c( "string"))


######################################
# TABLE 6
######################################

# Balance statistics
bal.pre <- r3(summary(m.out,standardize=T)$sum.all)[,c("Means Treated","Means Control","Std. Mean Diff.")]
bal.post <- r3(summary(m.out,standardize=T)$sum.matched)[,c("Means Treated","Means Control","Std. Mean Diff.")]
labz <- c("Propensity score","People killed (pre-T)","Partisan coercion (pre-T)","Partisan interdiction (pre-T)","People killed (3 months pre-T)","Partisan coercion (3 months pre-T)","Partisan interdiction (3 months pre-T)","Rail network centrality","Distance to rail station","Urbanization","Open terrain","Percent Jewish","Partisan control","Distance from frontline","District ID")
bal.mat <- rbind(rep("",3),bal.pre,rep("",3),bal.post)
bal.mat <- cbind(data.frame(V=c(paste0("Full dataset. ","N=",sum(summary(m.out)$nn[1,])/2," (T: ",summary(m.out)$nn[1,2]/2,", C: ",summary(m.out)$nn[1,1]/2,")"),labz,paste0("Matched dataset. ","N=",sum(summary(m.out)$nn[2,])/2," (T: ",summary(m.out)$nn[2,2]/2,", C: ",summary(m.out)$nn[2,1]/2,")"),labz)),bal.mat)
print.xtable(xtable(bal.mat),file="Analysis/Output2/Table_balance.tex",include.rownames=FALSE)


## Sensitivity Analysis: alternate time windows
## Rayon-level
load("Data/BELARUS_RAYONS_MONTHS_Full.RData")
load("Data/BELARUS_RAYONS_Full_sp.RData")
map_crd <- data.frame(RAYON_ID=map$RAYON_ID,LONG=coordinates(map)[,1],LAT=coordinates(map)[,2])
rayons.panel <- merge(rayons.panel,map_crd,by="RAYON_ID")
rayons.panel <- merge(rayons.panel,data.frame(YRMO=sort(unique(rayons.panel$YRMO)),TID=1:length(unique(rayons.panel$YRMO))),by="YRMO")
X <- rayons.panel
head(X)
drange <- 194308:194310
disag <- unique(map$RAYON_ID)

# Pre/post variables
i <- 23
dim(X)
mat.list <- lapply(1:length(disag),function(i){#print(paste(tvar,i))
  sub. <- X[X$RAYON_ID==disag[i],]
  tidz <- unique(sub.$TID[sub.$YRMO%in%drange])
  sub. <- sub.[order(sub.$TID),]
  sub. <- sub.[!duplicated(sub.$TID,fromLast=T),]
  sub2 <- data.frame(RAYON_ID=disag[i])
  sub2$P_PRE <- sum(sub.$PEOPLE_KILLED[sub.$TID<min(tidz)])
  sub2$P_PRE1 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$P_PRE2 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$P_PRE3 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$P_PRE4 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$P_PRE5 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$P_PRE6 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$P_DURING <- sum(sub.$PEOPLE_KILLED[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$P_POST <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)])
  sub2$P_POST1 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+1)])
  sub2$P_POST2 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+2)])
  sub2$P_POST3 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+3)])
  sub2$P_POST4 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+4)])
  sub2$P_POST5 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+5)])
  sub2$P_POST6 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+6)])
  sub2$H_PRE <- sum(sub.$HOUSES_DESTROYED[sub.$TID<min(tidz)])
  sub2$H_PRE1 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$H_PRE2 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$H_PRE3 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$H_PRE4 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$H_PRE5 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$H_PRE6 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$H_DURING <- sum(sub.$HOUSES_DESTROYED[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$H_POST <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)])
  sub2$H_POST1 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+1)])
  sub2$H_POST2 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+2)])
  sub2$H_POST3 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+3)])
  sub2$H_POST4 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+4)])
  sub2$H_POST5 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+5)])
  sub2$H_POST6 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+6)])
  sub2$R_PRE <- sum(sub.$REPRISAL[sub.$TID<min(tidz)])
  sub2$R_PRE1 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$R_PRE2 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$R_PRE3 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$R_PRE4 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$R_PRE5 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$R_PRE6 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$R_DURING <- sum(sub.$REPRISAL[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$R_POST <- sum(sub.$REPRISAL[sub.$TID>max(tidz)])
  sub2$R_POST1 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+1)])
  sub2$R_POST2 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+2)])
  sub2$R_POST3 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+3)])
  sub2$R_POST4 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+4)])
  sub2$R_POST5 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+5)])
  sub2$R_POST6 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+6)])
  sub2$T_CONVOYS_PRE <- sum(sub.$CONVOYS[sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE1 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE2 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE3 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE4 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE5 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE6 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_DURING <- sum(sub.$CONVOYS[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$T_COERCION_PRE <- sum(sub.$COERCION[sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE1 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE2 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE3 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE4 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE5 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE6 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$T_COERCION_DURING <- sum(sub.$COERCION[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$DISTLIB_MIN_PRE <- mean(sub.$DISTLIB_MIN[sub.$TID==min(tidz)])
  sub2
})
mat <- do.call(rbind,mat.list)
head(mat)

## Diff-in-diff
mat$T_CONVOYS <- 1*(mat$T_CONVOYS_DURING>mean(mat$T_CONVOYS_DURING))
mat$T_COERCION <- 1*(mat$T_COERCION_DURING>mean(mat$T_COERCION_DURING))
mat$ACTIVE <- 1*(mat$T_COERCION+mat$T_CONVOYS>0)
mat_pre <- mat
mat_pre$POST <- 0
mat_pre$P_ALL <- mat_pre$P_PRE
mat_pre$P_ALL1 <- mat_pre$P_PRE1
mat_pre$P_ALL2 <- mat_pre$P_PRE2
mat_pre$P_ALL3 <- mat_pre$P_PRE3
mat_pre$P_ALL4 <- mat_pre$P_PRE4
mat_pre$P_ALL5 <- mat_pre$P_PRE5
mat_pre$P_ALL6 <- mat_pre$P_PRE6
mat_pre$H_ALL <- mat_pre$H_PRE
mat_pre$H_ALL1 <- mat_pre$H_PRE1
mat_pre$H_ALL2 <- mat_pre$H_PRE2
mat_pre$H_ALL3 <- mat_pre$H_PRE3
mat_pre$H_ALL4 <- mat_pre$H_PRE4
mat_pre$H_ALL5 <- mat_pre$H_PRE5
mat_pre$H_ALL6 <- mat_pre$H_PRE6
mat_pre$R_ALL <- mat_pre$R_PRE
mat_pre$R_ALL1 <- mat_pre$R_PRE1
mat_pre$R_ALL2 <- mat_pre$R_PRE2
mat_pre$R_ALL3 <- mat_pre$R_PRE3
mat_pre$R_ALL4 <- mat_pre$R_PRE4
mat_pre$R_ALL5 <- mat_pre$R_PRE5
mat_pre$R_ALL6 <- mat_pre$R_PRE6
mat_pre$P_DALL <- mat_pre$P_PRE
mat_pre$P_DALL1 <- mat_pre$P_PRE1
mat_pre$P_DALL2 <- mat_pre$P_PRE2
mat_pre$P_DALL3 <- mat_pre$P_PRE3
mat_pre$P_DALL4 <- mat_pre$P_PRE4
mat_pre$P_DALL5 <- mat_pre$P_PRE5
mat_pre$P_DALL6 <- mat_pre$P_PRE6
mat_pre$H_DALL <- mat_pre$H_PRE
mat_pre$H_DALL1 <- mat_pre$H_PRE1
mat_pre$H_DALL2 <- mat_pre$H_PRE2
mat_pre$H_DALL3 <- mat_pre$H_PRE3
mat_pre$H_DALL4 <- mat_pre$H_PRE4
mat_pre$H_DALL5 <- mat_pre$H_PRE5
mat_pre$H_DALL6 <- mat_pre$H_PRE6
mat_pre$R_DALL <- mat_pre$R_PRE
mat_pre$R_DALL1 <- mat_pre$R_PRE1
mat_pre$R_DALL2 <- mat_pre$R_PRE2
mat_pre$R_DALL3 <- mat_pre$R_PRE3
mat_pre$R_DALL4 <- mat_pre$R_PRE4
mat_pre$R_DALL5 <- mat_pre$R_PRE5
mat_pre$R_DALL6 <- mat_pre$R_PRE6
mat_post <- mat
mat_post$POST <- 1
mat_post$P_ALL <- mat_post$P_POST
mat_post$P_ALL1 <- mat_post$P_POST1
mat_post$P_ALL2 <- mat_post$P_POST2
mat_post$P_ALL3 <- mat_post$P_POST3
mat_post$P_ALL4 <- mat_post$P_POST4
mat_post$P_ALL5 <- mat_post$P_POST5
mat_post$P_ALL6 <- mat_post$P_POST6
mat_post$H_ALL <- mat_post$H_POST
mat_post$H_ALL1 <- mat_post$H_POST1
mat_post$H_ALL2 <- mat_post$H_POST2
mat_post$H_ALL3 <- mat_post$H_POST3
mat_post$H_ALL4 <- mat_post$H_POST4
mat_post$H_ALL5 <- mat_post$H_POST5
mat_post$H_ALL6 <- mat_post$H_POST6
mat_post$R_ALL <- mat_post$R_POST
mat_post$R_ALL1 <- mat_post$R_POST1
mat_post$R_ALL2 <- mat_post$R_POST2
mat_post$R_ALL3 <- mat_post$R_POST3
mat_post$R_ALL4 <- mat_post$R_POST4
mat_post$R_ALL5 <- mat_post$R_POST5
mat_post$R_ALL6 <- mat_post$R_POST6
mat_post$P_DALL <- mat_post$P_POST+mat_post$P_DURING
mat_post$P_DALL1 <- mat_post$P_POST1+mat_post$P_DURING
mat_post$P_DALL2 <- mat_post$P_POST2+mat_post$P_DURING
mat_post$P_DALL3 <- mat_post$P_POST3+mat_post$P_DURING
mat_post$P_DALL4 <- mat_post$P_POST4+mat_post$P_DURING
mat_post$P_DALL5 <- mat_post$P_POST5+mat_post$P_DURING
mat_post$P_DALL6 <- mat_post$P_POST6+mat_post$P_DURING
mat_post$H_DALL <- mat_post$H_POST+mat_post$P_DURING
mat_post$H_DALL1 <- mat_post$H_POST1+mat_post$P_DURING
mat_post$H_DALL2 <- mat_post$H_POST2+mat_post$P_DURING
mat_post$H_DALL3 <- mat_post$H_POST3+mat_post$P_DURING
mat_post$H_DALL4 <- mat_post$H_POST4+mat_post$P_DURING
mat_post$H_DALL5 <- mat_post$H_POST5+mat_post$P_DURING
mat_post$H_DALL6 <- mat_post$H_POST6+mat_post$P_DURING
mat_post$R_DALL <- mat_post$R_POST+mat_post$P_DURING
mat_post$R_DALL1 <- mat_post$R_POST1+mat_post$P_DURING
mat_post$R_DALL2 <- mat_post$R_POST2+mat_post$P_DURING
mat_post$R_DALL3 <- mat_post$R_POST3+mat_post$P_DURING
mat_post$R_DALL4 <- mat_post$R_POST4+mat_post$P_DURING
mat_post$R_DALL5 <- mat_post$R_POST5+mat_post$P_DURING
mat_post$R_DALL6 <- mat_post$R_POST6+mat_post$P_DURING
mat_dind <- rbind(mat_pre,mat_post)
mat_dind <- mat_dind[mat_dind$ACTIVE==1,]


## Village-level
load("Data/BELARUS_VILLAGES_MONTHS_Full.RData")
load("Data/BELARUS_VILLAGES_Full.RData")
head(villages.panel)
villages.panel <- merge(villages.panel,data.frame(YRMO=sort(unique(villages.panel$YRMO)),TID=1:length(unique(villages.panel$YRMO))),by="YRMO")
X <- villages.panel
drange <- 194308:194310
disag <- unique(villages$VID)

# Pre/post variables
i <- 23
dim(X)
mat.list <- lapply(1:length(disag),function(i){print(i)
  sub. <- X[X$VID==disag[i],]
  tidz <- unique(sub.$TID[sub.$YRMO%in%drange])
  sub. <- sub.[order(sub.$TID),]
  sub. <- sub.[!duplicated(sub.$TID,fromLast=T),]
  sub2 <- data.frame(VID=disag[i])
  sub2$P_PRE <- sum(sub.$PEOPLE_KILLED[sub.$TID<min(tidz)])
  sub2$P_PRE1 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$P_PRE2 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$P_PRE3 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$P_PRE4 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$P_PRE5 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$P_PRE6 <- sum(sub.$PEOPLE_KILLED[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$P_DURING <- sum(sub.$PEOPLE_KILLED[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$P_POST <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)])
  sub2$P_POST1 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+1)])
  sub2$P_POST2 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+2)])
  sub2$P_POST3 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+3)])
  sub2$P_POST4 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+4)])
  sub2$P_POST5 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+5)])
  sub2$P_POST6 <- sum(sub.$PEOPLE_KILLED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+6)])
  sub2$H_PRE <- sum(sub.$HOUSES_DESTROYED[sub.$TID<min(tidz)])
  sub2$H_PRE1 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$H_PRE2 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$H_PRE3 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$H_PRE4 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$H_PRE5 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$H_PRE6 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$H_DURING <- sum(sub.$HOUSES_DESTROYED[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$H_POST <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)])
  sub2$H_POST1 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+1)])
  sub2$H_POST2 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+2)])
  sub2$H_POST3 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+3)])
  sub2$H_POST4 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+4)])
  sub2$H_POST5 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+5)])
  sub2$H_POST6 <- sum(sub.$HOUSES_DESTROYED[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+6)])
  sub2$R_PRE <- sum(sub.$REPRISAL[sub.$TID<min(tidz)])
  sub2$R_PRE1 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$R_PRE2 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$R_PRE3 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$R_PRE4 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$R_PRE5 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$R_PRE6 <- sum(sub.$REPRISAL[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$R_DURING <- sum(sub.$REPRISAL[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$R_POST <- sum(sub.$REPRISAL[sub.$TID>max(tidz)])
  sub2$R_POST1 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+1)])
  sub2$R_POST2 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+2)])
  sub2$R_POST3 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+3)])
  sub2$R_POST4 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+4)])
  sub2$R_POST5 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+5)])
  sub2$R_POST6 <- sum(sub.$REPRISAL[sub.$TID>max(tidz)&sub.$TID<=(max(tidz)+6)])
  sub2$T_CONVOYS_PRE <- sum(sub.$CONVOYS[sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE1 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE2 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE3 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE4 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE5 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_PRE6 <- sum(sub.$CONVOYS[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$T_CONVOYS_DURING <- sum(sub.$CONVOYS[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$T_COERCION_PRE <- sum(sub.$COERCION[sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE1 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-2)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE2 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-3)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE3 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-4)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE4 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-5)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE5 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-6)&sub.$TID<min(tidz)])
  sub2$T_COERCION_PRE6 <- sum(sub.$COERCION[sub.$TID>(min(tidz)-7)&sub.$TID<min(tidz)])
  sub2$T_COERCION_DURING <- sum(sub.$COERCION[sub.$TID>=min(tidz)&sub.$TID<=max(tidz)])
  sub2$DISTLIB_MIN_PRE <- mean(sub.$DISTLIB_MIN[sub.$TID==min(tidz)])
  sub2
})
mat <- do.call(rbind,mat.list)
head(mat)
save(mat, file="Data/BELARUS_VILLAGES_dindmat.RData")

## Diff-in-diff
load("Data/BELARUS_VILLAGES_dindmat.RData")
mat$T_CONVOYS <- 1*(mat$T_CONVOYS_DURING>mean(mat$T_CONVOYS_DURING))
mat$T_COERCION <- 1*(mat$T_COERCION_DURING>mean(mat$T_COERCION_DURING))
mat$T_CONVOYS <- 1*(mat$T_CONVOYS_DURING>0)
mat$T_COERCION <- 1*(mat$T_COERCION_DURING>0)
mat$ACTIVE <- 1*(mat$T_COERCION+mat$T_CONVOYS>0)
mat_pre <- mat
mat_pre$POST <- 0
mat_pre$P_ALL <- mat_pre$P_PRE
mat_pre$P_ALL1 <- mat_pre$P_PRE1
mat_pre$P_ALL2 <- mat_pre$P_PRE2
mat_pre$P_ALL3 <- mat_pre$P_PRE3
mat_pre$P_ALL4 <- mat_pre$P_PRE4
mat_pre$P_ALL5 <- mat_pre$P_PRE5
mat_pre$P_ALL6 <- mat_pre$P_PRE6
mat_pre$H_ALL <- mat_pre$H_PRE
mat_pre$H_ALL1 <- mat_pre$H_PRE1
mat_pre$H_ALL2 <- mat_pre$H_PRE2
mat_pre$H_ALL3 <- mat_pre$H_PRE3
mat_pre$H_ALL4 <- mat_pre$H_PRE4
mat_pre$H_ALL5 <- mat_pre$H_PRE5
mat_pre$H_ALL6 <- mat_pre$H_PRE6
mat_pre$R_ALL <- mat_pre$R_PRE
mat_pre$R_ALL1 <- mat_pre$R_PRE1
mat_pre$R_ALL2 <- mat_pre$R_PRE2
mat_pre$R_ALL3 <- mat_pre$R_PRE3
mat_pre$R_ALL4 <- mat_pre$R_PRE4
mat_pre$R_ALL5 <- mat_pre$R_PRE5
mat_pre$R_ALL6 <- mat_pre$R_PRE6
mat_pre$P_DALL <- mat_pre$P_PRE
mat_pre$P_DALL1 <- mat_pre$P_PRE1
mat_pre$P_DALL2 <- mat_pre$P_PRE2
mat_pre$P_DALL3 <- mat_pre$P_PRE3
mat_pre$P_DALL4 <- mat_pre$P_PRE4
mat_pre$P_DALL5 <- mat_pre$P_PRE5
mat_pre$P_DALL6 <- mat_pre$P_PRE6
mat_pre$H_DALL <- mat_pre$H_PRE
mat_pre$H_DALL1 <- mat_pre$H_PRE1
mat_pre$H_DALL2 <- mat_pre$H_PRE2
mat_pre$H_DALL3 <- mat_pre$H_PRE3
mat_pre$H_DALL4 <- mat_pre$H_PRE4
mat_pre$H_DALL5 <- mat_pre$H_PRE5
mat_pre$H_DALL6 <- mat_pre$H_PRE6
mat_pre$R_DALL <- mat_pre$R_PRE
mat_pre$R_DALL1 <- mat_pre$R_PRE1
mat_pre$R_DALL2 <- mat_pre$R_PRE2
mat_pre$R_DALL3 <- mat_pre$R_PRE3
mat_pre$R_DALL4 <- mat_pre$R_PRE4
mat_pre$R_DALL5 <- mat_pre$R_PRE5
mat_pre$R_DALL6 <- mat_pre$R_PRE6
mat_post <- mat
mat_post$POST <- 1
mat_post$P_ALL <- mat_post$P_POST
mat_post$P_ALL1 <- mat_post$P_POST1
mat_post$P_ALL2 <- mat_post$P_POST2
mat_post$P_ALL3 <- mat_post$P_POST3
mat_post$P_ALL4 <- mat_post$P_POST4
mat_post$P_ALL5 <- mat_post$P_POST5
mat_post$P_ALL6 <- mat_post$P_POST6
mat_post$H_ALL <- mat_post$H_POST
mat_post$H_ALL1 <- mat_post$H_POST1
mat_post$H_ALL2 <- mat_post$H_POST2
mat_post$H_ALL3 <- mat_post$H_POST3
mat_post$H_ALL4 <- mat_post$H_POST4
mat_post$H_ALL5 <- mat_post$H_POST5
mat_post$H_ALL6 <- mat_post$H_POST6
mat_post$R_ALL <- mat_post$R_POST
mat_post$R_ALL1 <- mat_post$R_POST1
mat_post$R_ALL2 <- mat_post$R_POST2
mat_post$R_ALL3 <- mat_post$R_POST3
mat_post$R_ALL4 <- mat_post$R_POST4
mat_post$R_ALL5 <- mat_post$R_POST5
mat_post$R_ALL6 <- mat_post$R_POST6
mat_post$P_DALL <- mat_post$P_POST+mat_post$P_DURING
mat_post$P_DALL1 <- mat_post$P_POST1+mat_post$P_DURING
mat_post$P_DALL2 <- mat_post$P_POST2+mat_post$P_DURING
mat_post$P_DALL3 <- mat_post$P_POST3+mat_post$P_DURING
mat_post$P_DALL4 <- mat_post$P_POST4+mat_post$P_DURING
mat_post$P_DALL5 <- mat_post$P_POST5+mat_post$P_DURING
mat_post$P_DALL6 <- mat_post$P_POST6+mat_post$P_DURING
mat_post$H_DALL <- mat_post$H_POST+mat_post$P_DURING
mat_post$H_DALL1 <- mat_post$H_POST1+mat_post$P_DURING
mat_post$H_DALL2 <- mat_post$H_POST2+mat_post$P_DURING
mat_post$H_DALL3 <- mat_post$H_POST3+mat_post$P_DURING
mat_post$H_DALL4 <- mat_post$H_POST4+mat_post$P_DURING
mat_post$H_DALL5 <- mat_post$H_POST5+mat_post$P_DURING
mat_post$H_DALL6 <- mat_post$H_POST6+mat_post$P_DURING
mat_post$R_DALL <- mat_post$R_POST+mat_post$P_DURING
mat_post$R_DALL1 <- mat_post$R_POST1+mat_post$P_DURING
mat_post$R_DALL2 <- mat_post$R_POST2+mat_post$P_DURING
mat_post$R_DALL3 <- mat_post$R_POST3+mat_post$P_DURING
mat_post$R_DALL4 <- mat_post$R_POST4+mat_post$P_DURING
mat_post$R_DALL5 <- mat_post$R_POST5+mat_post$P_DURING
mat_post$R_DALL6 <- mat_post$R_POST6+mat_post$P_DURING
mat_dind <- rbind(mat_pre,mat_post)
names(mat_dind)
mat_dind2 <- merge(mat_dind,villages,by="VID",all.x=T,all.y=F)
save(mat_dind2, file="Data/BELARUS_VILLAGES_dind.RData")
write.dta(mat_dind2, file="Data/BELARUS_VILLAGES_dind.dta",convert.factors = c( "string"))


## Matching
load("Data/BELARUS_VILLAGES_dind.RData")
tz <- 1:6
ti <- 2
for(ti in 1:6){print(ti)
  X <- mat_dind2[,c("VID","RAYON_ID","POST","T_CONVOYS","T_COERCION","P_ALL","H_ALL","R_ALL",paste0("P_ALL",tz[ti]),paste0("H_ALL",tz[ti]),paste0("R_ALL",tz[ti]),"P_DALL","H_DALL","R_DALL",paste0("P_DALL",tz[ti]),paste0("H_DALL",tz[ti]),paste0("R_DALL",tz[ti]),"P_DURING","H_DURING","R_DURING","P_PRE","H_PRE","R_PRE","T_CONVOYS_PRE","T_COERCION_PRE",paste0("P_PRE",tz[ti]),paste0("H_PRE",tz[ti]),paste0("R_PRE",tz[ti]),paste0("T_CONVOYS_PRE",tz[ti]),paste0("T_COERCION_PRE",tz[ti]),"CENT_BTWN_R","DIST_RAIL","URBAN100","OPEN_TERRAIN","ETH_JUD","CONTROL_AREA","DISTLIB_MIN_PRE")]
  
  # CONVOYS
  formz <- as.formula(paste0("T_CONVOYS ~ P_PRE + T_COERCION_PRE+ T_CONVOYS_PRE +P_PRE",tz[ti],"+ T_COERCION_PRE",tz[ti]," + T_CONVOYS_PRE",tz[ti],"  + CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN_PRE"))
  m.out <- matchit(formz, distance="logit",exact=c("RAYON_ID"), ratio=1, data=X, discard="both",reestimate=T,caliper=.25)
  m.data <- match.data(m.out)
  write.dta(m.data, file=paste0("Data/BELARUS_VILLAGES_matched_v1_tz",tz[ti],".dta"),convert.factors = c( "string"))
  save(m.data, file=paste0("Data/BELARUS_VILLAGES_matched_v1_tz",tz[ti],".RData"))
  
  # COERCION
  formz <- as.formula(paste0("T_COERCION ~ P_PRE + T_COERCION_PRE+ T_CONVOYS_PRE +P_PRE",tz[ti],"+ T_COERCION_PRE",tz[ti]," + T_CONVOYS_PRE",tz[ti],"  + CENT_BTWN_R+DIST_RAIL+URBAN100+OPEN_TERRAIN+ETH_JUD+CONTROL_AREA+DISTLIB_MIN_PRE"))
  m.out <- matchit(formz, distance="logit",exact=c("RAYON_ID"), ratio=1, data=X, discard="both",reestimate=T,caliper=.25)
  m.data <- match.data(m.out)
  write.dta(m.data, file=paste0("Data/BELARUS_VILLAGES_matched_v2_tz",tz[ti],".dta"),convert.factors = c( "string"))
  save(m.data, file=paste0("Data/BELARUS_VILLAGES_matched_v2_tz",tz[ti],".RData"))
}


#####################################
## to implement D-in-D,
## run repcode_table5.do in STATA
#####################################





#####################################
## FIGURE 6
#####################################

load("Data/BELARUS_RAYONS_MONTHS_Full.RData")
load("Data/BELARUS_RAYONS_Full_sp.RData")
map_crd <- data.frame(RAYON_ID=map$RAYON_ID,LONG=coordinates(map)[,1],LAT=coordinates(map)[,2])
rayons.panel <- merge(rayons.panel,map_crd,by="RAYON_ID")

## Weights + MSM
outz <- c("PEOPLE_KILLED","HOUSES_DESTROYED","REPRISAL","PEOPLE_KILLED_p1","HOUSES_DESTROYED_p1","REPRISAL_p1")
varz <- c("CONVOYS","COERCION","COERCION_CONVOYS")
k <- 1
o <- 4

## Iterate over outcomes
out.list <- lapply(1:length(outz),function(o){
  ## Iterate over treatment regimes
  te.list <- lapply(1:3,function(k){
    X <- rayons.panel
    X$TREAT_t1 <- streg.tlag(data=X, timevar="YRMO", lagvar=paste0("TREAT_",varz[k]), order=1)
    X$TREAT_t1[is.na(X$TREAT_t1)] <- 0
    print(paste(outz[o],varz[k]))
    # Weights
    f.n <- as.formula(paste0("TREAT_",varz[k]," ~ T_HIST_",varz[k]," + as.factor(RAYON_ID)"))
    f.d <- as.formula(paste0("TREAT_",varz[k]," ~ T_HIST_",varz[k]," + DISTLIB_MIN + as.factor(RAYON_ID)"))
    ps.n <- glm(f.n, data= X, family="binomial")
    ps.d <- glm(f.d, data= X, family="binomial")
    id.n <- row.names(ps.n$model)
    id.d <- row.names(ps.d$model)
    idz <- intersect(id.n,id.d)
    X <- X[row.names(X)%in%idz,]
    X$ps.n <- predict(ps.n,type="response")[row.names(ps.n$model)%in%idz]
    X$ps.d <- predict(ps.d,type="response")[row.names(ps.d$model)%in%idz]
    disag <- sort(unique(X$YRMO))
    ipw.list <- lapply(1:length(disag),function(t){
      sub.X <- X[X$YRMO<=disag[t],]  
      sub.X <- aggregate(sub.X[,c("ps.n","ps.d")],by=list(RAYON_ID=sub.X$RAYON_ID),prod)
      sub.X$YRMO <- disag[t]
      sub.X$IPW <- sub.X$ps.n/sub.X$ps.d
      #sub.X$IPW <- (sub.X$ps.n+1)/(sub.X$ps.d+1)
      sub.X <- sub.X[,c("RAYON_ID","YRMO","IPW")]
      sub.X$IPW[is.na(sub.X$IPW)|!is.finite(sub.X$IPW)] <- 0
      #print(paste(disag[t],max(sub.X$IPW)))
      names(sub.X)[names(sub.X)=="IPW"] <- paste0("IPW_",varz[k])
      sub.X
    })
    ipw.mat <- do.call(rbind,ipw.list)
    summary(ipw.mat)
    X <- merge(X,ipw.mat,by=c("RAYON_ID","YRMO"),all.x=F,all.y=T)
    X <- X[X$YRMO < 194407,]
    X <- merge(X,data.frame(YRMO=sort(unique(X$YRMO)),TID=1:length(unique(X$YRMO))),by="YRMO")
    
    ## ATE matrix
    temat <- data.frame(Estimator=c("Contemporaneous treatment effect","Treatment history effect","All or nothing effect","Local treatment assignment effect","Neighbor effect on controls","Neighbor effect on treated","Interaction"),Estimate=NA,Lower=NA,Upper=NA)
    ## MSM estimator
    fmz <- as.formula(paste0("I(log(",outz[o],"+1)) ~ as.factor(YRMO) + as.factor(RAYON_ID) + s(TID) + TREAT_",varz[k]," + T_HIST_",varz[k]))
    w <- X[,paste0("IPW_",varz[k])]
    mod <- gam(fmz, weights=w, data=X);summary(mod)
    betahat <- mod$coef
    betahat <- betahat[!is.na(betahat)]
    vcmat <- summary(mod)$cov.scaled
    vcmat <- vcmat[row.names(vcmat)%in%names(betahat),colnames(vcmat)%in%names(betahat)]
    betas <- rmvnorm(n=1000, mean = betahat, sigma = vcmat, method="svd")
    b_treat <- betas[,paste0("TREAT_",varz[k])]
    b_thist <- betas[,paste0("T_HIST_",varz[k])]
    
    # Contemporaneous treatment effect
    ezt <- b_treat
    temat[1,c("Estimate","Lower","Upper")] <- c(mean(ezt),quantile(ezt,.025),quantile(ezt,.975))
    
    # Treatment history effect
    ezt <- b_thist
    temat[2,c("Estimate","Lower","Upper")] <- c(mean(ezt),quantile(ezt,.025),quantile(ezt,.975))
    
    
    ## Spatial estimator
    X$INTX <- X[,paste0("TREAT_",varz[k])]*X[,paste0("W_TREAT_",varz[k])]
    fmz <- as.formula(paste0("I(log(",outz[o],"+1)) ~ as.factor(YRMO) +as.factor(RAYON_ID) + TREAT_",varz[k]," + T_HIST_",varz[k]," + W_TREAT_",varz[k],"+ INTX"))
    w <- X[,paste0("IPW_",varz[k])]
    mod <- gam(fmz, weights=w, data=X);summary(mod)
    betahat <- mod$coef
    betahat <- betahat[!is.na(betahat)]
    vcmat <- summary(mod)$cov.scaled
    vcmat <- vcmat[row.names(vcmat)%in%names(betahat),colnames(vcmat)%in%names(betahat)]
    betas <- rmvnorm(n=1000, mean = betahat, sigma = vcmat, method="svd")
    b_treat <- betas[,paste0("TREAT_",varz[k])]
    b_wt <- betas[,paste0("W_TREAT_",varz[k])]
    b_twt <- 0
    if(length(which(colnames(betas)%in%"INTX"))>0){b_twt <- betas[,"INTX"]}
    b_thist <- betas[,paste0("T_HIST_",varz[k])]
    
    # All or nothing
    ezt <- b_treat + 1*(b_wt+b_twt)
    temat[3,c("Estimate","Lower","Upper")] <- c(mean(ezt),quantile(ezt,.025),quantile(ezt,.975))
    
    # Local treatment assignment
    ezt <- b_treat + .5*(b_twt)
    temat[4,c("Estimate","Lower","Upper")] <- c(mean(ezt),quantile(ezt,.025),quantile(ezt,.975))
    
    # Neighbor effect on controls
    ezt <- .5*(b_wt)
    temat[5,c("Estimate","Lower","Upper")] <- c(mean(ezt),quantile(ezt,.025),quantile(ezt,.975))
    
    # Neighbor effect on treated
    ezt <- .5*(b_wt+b_twt)
    temat[6,c("Estimate","Lower","Upper")] <- c(mean(ezt),quantile(ezt,.025),quantile(ezt,.975))
    
    # Diff in neighbor effect on treated vs. control
    ezt <- .5*(b_twt)
    temat[7,c("Estimate","Lower","Upper")] <- c(mean(ezt),quantile(ezt,.025),quantile(ezt,.975))
    
    temat
  })
  names(te.list) <- varz
  te.list
})
names(out.list) <- outz
out.list



## FIGURE 6a
# Contemporaneous
xmin <- -65
xmax <- 375
png(file="Output/FIGURE6a.png",width=8,height=6,units="in",res=300)
par(mar=c(2,6,0,0))
plot(x=0,y=0,col=NA,xlim=c(xmin,xmax),ylim=c(0,12),bty="n",xaxt="n",yaxt="n",ylab="",xlab="")
rect(xleft=xmin,xright=xmax,ybottom=10,ytop=12,col="gray80",border=NA)
rect(xleft=xmin,xright=xmax,ybottom=8,ytop=10,col="gray90",border=NA)
rect(xleft=xmin,xright=xmax,ybottom=6,ytop=8,col="gray80",border=NA)
rect(xleft=xmin,xright=xmax,ybottom=4,ytop=6,col="gray90",border=NA)
rect(xleft=xmin,xright=xmax,ybottom=2,ytop=4,col="gray80",border=NA)
rect(xleft=xmin,xright=xmax,ybottom=0,ytop=2,col="gray90",border=NA)
segments(x0=xmin,x1=xmax,y0=12,y1=12)
segments(x0=xmin,x1=xmax,y0=8,y1=8)
segments(x0=xmin,x1=xmax,y0=4,y1=4)
segments(x0=xmin,x1=xmax,y0=0,y1=0)
axis(1,at=pretty(xmin:xmax,n=10),cex.axis=1)
rect(xleft=100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS$Lower[1])-1),xright=100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS$Upper[1])-1),ybottom=10.5,ytop=11.5,border=NA,col="purple")
segments(x0=100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS$Estimate[1])-1),x1=100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS$Estimate[1])-1),y0=10.4,y1=11.6,col="black",lwd=2)
rect(xleft=100*(exp(out.list$PEOPLE_KILLED_p1$COERCION$Lower[1])-1),xright=100*(exp(out.list$PEOPLE_KILLED_p1$COERCION$Upper[1])-1),ybottom=8.5,ytop=9.5,border=NA,col="black")
segments(x0=100*(exp(out.list$PEOPLE_KILLED_p1$COERCION$Estimate[1])-1),x1=100*(exp(out.list$PEOPLE_KILLED_p1$COERCION$Estimate[1])-1),y0=8.4,y1=9.6,col="white",lwd=2)
rect(xleft=100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS$Lower[1])-1),xright=100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS$Upper[1])-1),ybottom=6.5,ytop=7.5,border=NA,col="purple")
segments(x0=100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS$Estimate[1])-1),x1=100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS$Estimate[1])-1),y0=6.4,y1=7.6,col="black",lwd=2)
rect(xleft=100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION$Lower[1])-1),xright=100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION$Upper[1])-1),ybottom=4.5,ytop=5.5,border=NA,col="black")
segments(x0=100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION$Estimate[1])-1),x1=100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION$Estimate[1])-1),y0=4.4,y1=5.6,col="white",lwd=2)
rect(xleft=100*(exp(out.list$REPRISAL_p1$CONVOYS$Lower[1])-1),xright=100*(exp(out.list$REPRISAL_p1$CONVOYS$Upper[1])-1),ybottom=2.5,ytop=3.5,border=NA,col="purple")
segments(x0=100*(exp(out.list$REPRISAL_p1$CONVOYS$Estimate[1])-1),x1=100*(exp(out.list$REPRISAL_p1$CONVOYS$Estimate[1])-1),y0=2.4,y1=3.6,col="black",lwd=2)
rect(xleft=100*(exp(out.list$REPRISAL_p1$COERCION$Lower[1])-1),xright=100*(exp(out.list$REPRISAL_p1$COERCION$Upper[1])-1),ybottom=0.5,ytop=1.5,border=NA,col="black")
segments(x0=100*(exp(out.list$REPRISAL_p1$COERCION$Estimate[1])-1),x1=100*(exp(out.list$REPRISAL_p1$COERCION$Estimate[1])-1),y0=0.4,y1=1.6,col="white",lwd=2)
segments(x0=0,x1=0,y0=-1,y1=12,col="red")
par(xpd="n")
text(x=xmin,y=10,pos=2,labels = "People\n killed\n ",cex=1.5)
text(x=xmin,y=6,pos=2,labels = "Houses\n destroyed\n ",cex=1.5)
text(x=xmin,y=2,pos=2,labels = "Number of\n reprisals\n ",cex=1.5)
text(x=xmin,y=10,pos=2,labels = "\n\n (% change)",cex=1.2,col="darkgrey")
text(x=xmin,y=6,pos=2,labels = "\n\n (% change)",cex=1.2,col="darkgrey")
text(x=xmin,y=2,pos=2,labels = "\n\n (% change)",cex=1.2,col="darkgrey")
text(x=max(100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS$Upper[1])-1),0),y=11,"Interdiction\n(convoy derailment)",pos=4,col="purple4",cex=1.2)
text(x=100*(exp(out.list$PEOPLE_KILLED_p1$COERCION$Upper[1])-1),y=9,"Coercion\n(garrison attack)",pos=4,col="black",cex=1.2)
dev.off()

## FIGURE 6b
# Cumulative
xmin <- -6
xmax <- 24
png(file="Output/FIGURE6b.png",width=8,height=6,units="in",res=300)
par(mar=c(2,6,0,0))
plot(x=0,y=0,col=NA,xlim=c(xmin,xmax),ylim=c(0,12),bty="n",xaxt="n",yaxt="n",ylab="",xlab="")
rect(xleft=xmin,xright=xmax,ybottom=10,ytop=12,col="gray80",border=NA)
rect(xleft=xmin,xright=xmax,ybottom=8,ytop=10,col="gray90",border=NA)
rect(xleft=xmin,xright=xmax,ybottom=6,ytop=8,col="gray80",border=NA)
rect(xleft=xmin,xright=xmax,ybottom=4,ytop=6,col="gray90",border=NA)
rect(xleft=xmin,xright=xmax,ybottom=2,ytop=4,col="gray80",border=NA)
rect(xleft=xmin,xright=xmax,ybottom=0,ytop=2,col="gray90",border=NA)
segments(x0=xmin,x1=xmax,y0=12,y1=12)
segments(x0=xmin,x1=xmax,y0=8,y1=8)
segments(x0=xmin,x1=xmax,y0=4,y1=4)
segments(x0=xmin,x1=xmax,y0=0,y1=0)
axis(1,at=pretty(xmin:xmax,n=10),cex.axis=1)
rect(xleft=100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS$Lower[2])-1),xright=100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS$Upper[2])-1),ybottom=10.5,ytop=11.5,border=NA,col="purple")
segments(x0=100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS$Estimate[2])-1),x1=100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS$Estimate[2])-1),y0=10.4,y1=11.6,col="black",lwd=2)
rect(xleft=100*(exp(out.list$PEOPLE_KILLED_p1$COERCION$Lower[2])-1),xright=100*(exp(out.list$PEOPLE_KILLED_p1$COERCION$Upper[2])-1),ybottom=8.5,ytop=9.5,border=NA,col="black")
segments(x0=100*(exp(out.list$PEOPLE_KILLED_p1$COERCION$Estimate[2])-1),x1=100*(exp(out.list$PEOPLE_KILLED_p1$COERCION$Estimate[2])-1),y0=8.4,y1=9.6,col="white",lwd=2)
rect(xleft=100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS$Lower[2])-1),xright=100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS$Upper[2])-1),ybottom=6.5,ytop=7.5,border=NA,col="purple")
segments(x0=100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS$Estimate[2])-1),x1=100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS$Estimate[2])-1),y0=6.4,y1=7.6,col="black",lwd=2)
rect(xleft=100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION$Lower[2])-1),xright=100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION$Upper[2])-1),ybottom=4.5,ytop=5.5,border=NA,col="black")
segments(x0=100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION$Estimate[2])-1),x1=100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION$Estimate[2])-1),y0=4.4,y1=5.6,col="white",lwd=2)
rect(xleft=100*(exp(out.list$REPRISAL_p1$CONVOYS$Lower[2])-1),xright=100*(exp(out.list$REPRISAL_p1$CONVOYS$Upper[2])-1),ybottom=2.5,ytop=3.5,border=NA,col="purple")
segments(x0=100*(exp(out.list$REPRISAL_p1$CONVOYS$Estimate[2])-1),x1=100*(exp(out.list$REPRISAL_p1$CONVOYS$Estimate[2])-1),y0=2.4,y1=3.6,col="black",lwd=2)
rect(xleft=100*(exp(out.list$REPRISAL_p1$COERCION$Lower[2])-1),xright=100*(exp(out.list$REPRISAL_p1$COERCION$Upper[2])-1),ybottom=0.5,ytop=1.5,border=NA,col="black")
segments(x0=100*(exp(out.list$REPRISAL_p1$COERCION$Estimate[2])-1),x1=100*(exp(out.list$REPRISAL_p1$COERCION$Estimate[2])-1),y0=0.4,y1=1.6,col="white",lwd=2)
segments(x0=0,x1=0,y0=-1,y1=12,col="red")
par(xpd="n")
text(x=xmin,y=10,pos=2,labels = "People\n killed\n ",cex=1.5)
text(x=xmin,y=6,pos=2,labels = "Houses\n destroyed\n ",cex=1.5)
text(x=xmin,y=2,pos=2,labels = "Number of\n reprisals\n ",cex=1.5)
text(x=xmin,y=10,pos=2,labels = "\n\n (% change)",cex=1.2,col="darkgrey")
text(x=xmin,y=6,pos=2,labels = "\n\n (% change)",cex=1.2,col="darkgrey")
text(x=xmin,y=2,pos=2,labels = "\n\n (% change)",cex=1.2,col="darkgrey")
text(x=max(100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS$Upper[2])-1),0),y=11,"Interdiction\n(convoy derailment)",pos=4,col="purple4",cex=1.2)
text(x=100*(exp(out.list$PEOPLE_KILLED_p1$COERCION$Upper[2])-1),y=9,"Coercion\n(garrison attack)",pos=4,col="black",cex=1.2)
dev.off()


#########################################
## CFs
#########################################


# Contemporaneous treatment effect
100*(exp(out.list$PEOPLE_KILLED$CONVOYS[1,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION[1,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION_CONVOYS[1,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS[1,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION[1,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION_CONVOYS[1,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED$CONVOYS[1,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION[1,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION_CONVOYS[1,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS[1,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION[1,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION_CONVOYS[1,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL$CONVOYS[1,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION[1,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION_CONVOYS[1,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL_p1$CONVOYS[1,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION[1,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION_CONVOYS[1,c("Estimate","Lower","Upper")])-1)

# Treatment history effect
100*(exp(out.list$PEOPLE_KILLED$CONVOYS[2,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION[2,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION_CONVOYS[2,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS[2,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION[2,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION_CONVOYS[2,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED$CONVOYS[2,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION[2,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION_CONVOYS[2,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS[2,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION[2,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION_CONVOYS[2,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL$CONVOYS[2,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION[2,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION_CONVOYS[2,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL_p1$CONVOYS[2,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION[2,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION_CONVOYS[2,c("Estimate","Lower","Upper")])-1)


# All or nothing effect
100*(exp(out.list$PEOPLE_KILLED$CONVOYS[3,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION[3,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION_CONVOYS[3,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS[3,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION[3,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION_CONVOYS[3,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED$CONVOYS[3,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION[3,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION_CONVOYS[3,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS[3,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION[3,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION_CONVOYS[3,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL$CONVOYS[3,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION[3,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION_CONVOYS[3,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL_p1$CONVOYS[3,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION[3,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION_CONVOYS[3,c("Estimate","Lower","Upper")])-1)


# Local treatment assignment effect
100*(exp(out.list$PEOPLE_KILLED$CONVOYS[4,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION[4,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION_CONVOYS[4,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS[4,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION[4,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION_CONVOYS[4,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED$CONVOYS[4,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION[4,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION_CONVOYS[4,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS[4,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION[4,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION_CONVOYS[4,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL$CONVOYS[4,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION[4,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION_CONVOYS[4,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL_p1$CONVOYS[4,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION[4,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION_CONVOYS[4,c("Estimate","Lower","Upper")])-1)


# Neighbor effect on controls
100*(exp(out.list$PEOPLE_KILLED$CONVOYS[5,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION[5,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION_CONVOYS[5,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS[5,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION[5,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION_CONVOYS[5,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED$CONVOYS[5,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION[5,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION_CONVOYS[5,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS[5,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION[5,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION_CONVOYS[5,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL$CONVOYS[5,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION[5,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION_CONVOYS[5,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL_p1$CONVOYS[5,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION[5,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION_CONVOYS[5,c("Estimate","Lower","Upper")])-1)


# Neighbor effect on treated
100*(exp(out.list$PEOPLE_KILLED$CONVOYS[6,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION[6,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION_CONVOYS[6,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS[6,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION[6,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION_CONVOYS[6,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED$CONVOYS[6,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION[6,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION_CONVOYS[6,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS[6,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION[6,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION_CONVOYS[6,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL$CONVOYS[6,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION[6,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION_CONVOYS[6,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL_p1$CONVOYS[6,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION[6,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION_CONVOYS[6,c("Estimate","Lower","Upper")])-1)


# Interaction
100*(exp(out.list$PEOPLE_KILLED$CONVOYS[7,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION[7,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED$COERCION_CONVOYS[7,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$PEOPLE_KILLED_p1$CONVOYS[7,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION[7,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$PEOPLE_KILLED_p1$COERCION_CONVOYS[7,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED$CONVOYS[7,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION[7,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED$COERCION_CONVOYS[7,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$HOUSES_DESTROYED_p1$CONVOYS[7,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION[7,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$HOUSES_DESTROYED_p1$COERCION_CONVOYS[7,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL$CONVOYS[7,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION[7,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL$COERCION_CONVOYS[7,c("Estimate","Lower","Upper")])-1)

100*(exp(out.list$REPRISAL_p1$CONVOYS[7,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION[7,c("Estimate","Lower","Upper")])-1)
100*(exp(out.list$REPRISAL_p1$COERCION_CONVOYS[7,c("Estimate","Lower","Upper")])-1)
